Structure–composition trends in multicomponent borosilicate-based glasses deduced from molecular dynamics simulations with improved B–O and P–O force fields

Abstract

We present a comprehensive molecular dynamics (MD) simulation study of composition–structure trends in a set of 25 glasses of widely spanning compositions from the following four systems of increasing complexity: Na₂O–B₂O₃, Na₂O–B₂O₃–SiO₂, Na₂O–CaO–SiO₂–P₂O₅, and Na₂O–CaO–B₂O₃–SiO₂–P₂O₅. The simulations involved new B–O and P–O potential parameters developed within the polarizable shell-model framework, thereby combining the beneficial features of an overall high accuracy and excellent transferability among different glass systems and compositions: this was confirmed by the good accordance with experimental data on the relative BO₃/BO₄ populations in borate and boro(phospho)silicate networks, as well as with the orthophosphate fractions in bioactive (boro)phosphosilicate glasses, which is believed to strongly influence their bone-bonding properties. The bearing of the simulated melt-cooling rate on the borate/phosphate speciations is discussed. Each local {BO₃, BO₄, SiO₄, PO₄} coordination environment remained independent of the precise set of co-existing network formers, while all trends observed in bond-lengths/angles mainly reflected the glass-network polymerization, i.e., the relative amounts of bridging oxygen (BO) and non-bridging oxygen (NBO) species. The structural roles of the Na⁺/Ca²⁺ cations were also probed, targeting their local coordination environments and their relative preferences to associate with the various borate, silicate, and phosphate moieties. We evaluate and discuss the common classification of alkali/alkaline-earth metal ions as charge-compensators of either BO₄ tetrahedra or NBO anions in borosilicate glasses, also encompassing the less explored NBO-rich regime: the Na⁺/Ca²⁺ cations mainly associate with BO/NBO species of SiO₄/BO₃ groups, with significant relative Na–BO₄ contacts only observed in B-rich glass networks devoid of NBO species, whereas NBO-rich glass networks also reveal substantial amounts of NBO-bearing BO₄ tetrahedra.

---

1 Introduction

Underlying the exploitation of glass for numerous purposes—ranging from everyday products to advanced industrial applications—is the well-controlled tuning of the desirable physical and chemical properties often attainable by a straightforward variation of the oxide precursor contents. This tailoring is greatly assisted by a detailed insight into the underlying composition–structure relations dictating the targeted glass properties. Yet, the amorphous nature yields distributions of the experimental observables, which limits the experimental probing and thereby also the structural insight, notably so at the medium-range scale (0.3–1 nm). The limited detailed structural understanding hampers a rational glass design for establishing composition–property trends, usually necessitating empirical “trial-and-error” procedures.

Here we consider oxide-based glasses involving at least one of the network formers B, Si, and P, which encompass all primary formers but Al, thereby accounting for large groups of glasses of importance to chemistry, physics, as well as to materials and earth Sciences. The glass network-forming species are interlinked by covalent –O– bridges.^1,2 The oxygen speciation is dominated either by such bridging oxygen (BO) atoms or by non-bridging oxygen (NBO) ions, which coordinate two (O^[2]) and one (O^[1]) network formers, respectively. The most frequent cation coordination number is four, as commonly adopted by B (BO₄; B^[4] coordination) and exclusively by Si (SiO₄) and P (PO₄).^1,2 However, the small B atom often forms triangular BO₃ groups (B^[3] coordinations), as for instance encountered in vitreous B₂O₃.^1–3

Electropositive cations (e.g., Na⁺ and Ca²⁺) are often referred to as glass-network modifiers, because they depolymerize the networks by arranging BO → NBO conversions.^1,2 The modifiers mainly balance the negative charges of NBO ions and other anionic species, such as [BO₄]⁻ and phosphate groups present in B and P bearing glasses; notably, a commonly adopted (albeit simplified) structural description of borate/borosilicate glasses assumes that the M⁺/M²⁺ cations act predominantly either as NBO-associated “modifiers” or as “compensators” of the [BO₄]⁻ tetrahedra, as most frequently discussed for Na₂O–B₂O₃–SiO₂ glasses.^4–10 Herein, we examine the relevance of this rather categorical classification in various B-based glass systems, also targeting the less investigated scenario of NBO-rich networks.

Concerning the experimentally derived structural insight from glasses involving B, Si, and/or P as network formers, diffraction techniques offer average cation–oxygen distances and bond-angle distributions (BADs),^11–18 while spectroscopic methods—particularly magic-angle spinning (MAS) nuclear magnetic resonance (NMR)—provide accurate relative populations of the BO₃ and BO₄ coordinations (denoted by x^[3]_B and x^[4]_B, respectively) in B-bearing glasses.^{2,4–8,19–33} In phosphosilicate-based glasses, the distributions of co-existing {Qⁿ_P} and {Qⁿ_Si} tetrahedra with varying number n of BO atoms at the respective PO₄ and SiO₄ tetrahedra are often attainable by ³¹P/²⁹Si MAS NMR,^2,30,32–36 whereas ¹⁷O NMR informs about the BO/NBO speciation in oxide glasses.^{2,21,22,37,38} Solid-state NMR may also unveil the medium-range glass structure, such as providing connectivities/proximities and spatial distributions among the various network formers,^{2,28,31,39–44} but sophisticated and time-consuming experimentation is required that remains sparsely exploited for multicomponent glasses as compared with routine ¹¹B/²⁹Si/³¹P MAS NMR applications.

Given the limited accessible experimental data on several structural features of glasses, atomistic molecular dynamics (MD) simulations offer an alternative probing up to the nanometer scale,^{10,34,36,42–62} where a complete set of structural parameters is available from one sole glass model. Yet, the predictive power of classical MD simulations depends critically on the precise choice of force fields to model the local cation–oxygen interactions, which ultimately also control the medium-range organization of atoms. For instance, the temperature dependent BO₃ ⇄ BO₄ equilibrium in melts shifts to the right during melt-cooling, thereby elevating the BO₄ content in the glass structure relative to the melt;^{8–10,23,24,26} this feature coupled with a strong composition dependence of the borate speciation compromises the accuracy of modeled {BO₃, BO₄} populations in B-bearing glasses. These problems were mitigated by invoking glass-composition dependent B–O potential parameters,^52–54 however at the expense of system-dependent force fields with limited transferability.

Force fields accounting for polarization effects are known to enhance the modeling of both short and medium range structures of B₂O₃, encompassing boroxol-ring formation^46,63–65 and the presence of BO₄ groups at extreme pressures.⁶⁶ However, thus far only two such force-field options exist for B-based multicomponent glasses.^10,44 Here we present the derivation and assessments of new B–O and P–O pair-potentials utilized by us in a recent structural report on borophosphosilicate (BPS) glasses.⁴⁴ They were developed within the polarizable shell-model potential introduced by Sanders et al.⁴⁵ for the modeling of SiO₂, which was subsequently extended to Na₂O–CaO–SiO₂–P₂O₅ glasses.^55–57 The shell model employs full formal cation charges, which facilitates its implementation in widely spanning classes of multicomponent glasses. The high predictive abilities of our proposed force fields are demonstrated by validations on crystalline B/P-bearing structures (see the ESI†), as well as by the structural probing of glasses from the Na₂O–B₂O₃, Na₂O–B₂O₃–SiO₂, Na₂O–CaO–SiO₂–P₂O₅, and Na₂O–CaO–B₂O₃–SiO₂–P₂O₅ systems.

Phosphorus participates only to a limited extent in the Na₂O–CaO–(B₂O₃)–SiO₂–P₂O₅ glass networks targeted herein that belong to the group of bioactive glasses (BGs) with applications to bone grafting and tissue engineering.^67,68 Their bone-bonding properties correlate qualitatively with the extent of hydroxyapatite [Ca₅(PO₄)₃OH] mineralization observed in vitro, which is believed to mainly depend on the amount of ortho-phosphate (Q⁰_P) anions in the glass along with its silicate network connectivity;^69–73 the latter is denoted by N̄^Si_BO and represents the average number of BO atoms per SiO₄ tetrahedron.⁷³ Silicate-based BGs manifest comparatively fragmented glass networks conforming to 2.0 ≤ N̄^Si_BO ≤ 2.6,⁷³ while featuring sufficient amounts of Na⁺/Ca²⁺ cations to (potentially) charge-balance all phosphate species as Q⁰_P moieties. Notwithstanding that ³¹P MAS NMR confirms that Q⁰_P tetrahedra dominate the phosphate speciations (>85%), there are also non-negligible contributions from Q¹_P groups (signifying one P–O–Si bridge^{36,39–42,56,57}), whereas higher-polymerization Qⁿ_P (n ≥ 2) species are essentially absent.^36,41 Unfortunately, all force fields utilized thus far in classical MD simulations of BGs give relatively poor predictions of the {Qⁿ_P} populations, where the orthophosphate fraction is consistently underestimated relative to experiments,^34,36,42 even when employing the likely best P–O potential available.⁵⁷ Here we demonstrate that refinement of those parameters significantly improves the MD-derived {Qⁿ_P} speciations in phosphosilicate glass models throughout the composition range relevant for BGs.

This article is organized as follows: Section 2 introduces the glass systems/samples and common notation, while the MD simulation procedures and force-field developments are described in Section 3. Sections 4 and 5 contrast the MD-derived borate and phosphate speciations, respectively, against experimental results, including their dependence on the choice of melt-cooling rate in the simulation. Section 6 discusses the local coordination environments of the network formers F = {B^[3], B^[4], Si, P}, including the F–O distances and intra/interpolyhedral bond angles, as well as the relative propensities for F–NBO associations. Section 7 discusses the structural roles of the Na⁺/Ca²⁺ cations, encompassing their local coordination environments and their relative preferences to associate with the various NBO-bearing silicate, borate, and phosphate moieties, with a particular focus on quantifying the extents of “modifier/compensator” roles from a wider perspective than that normally adopted. Section 8 summarizes the main findings of our study.

2 Glass samples

Table 1 compiles our set of 25 modeled glasses, encompassing members from the Na₂O–B₂O₃ (NB), Na₂O–B₂O₃–SiO₂ (NBS), Na₂O–CaO–SiO₂–P₂O₅ (NCPS), and Na₂O–CaO–B₂O₃–SiO₂–P₂O₅ (NCBPS) systems, where each capital letter in the acronyms specifies the respective element of the glass-network modifier {Na, Ca} or former {B, P, Si}. We reserve the symbol “x” for fractions; for instance, the molar fraction of element E and oxide “EO” is denoted by x_E and x(EO), respectively. All NC(B)PS glass members considered herein were prepared and characterized in our previous work.^30,36,43,44 The borate/borosilicate glass compositions and their accompanying experimental data were obtained from the literature, as specified in Table 1 and Table S1 (ESI†).

Table 1 Glass samples and MD-derived data^a

Glass	Oxide equivalents (mol%)						BO4 fraction^c				Na/Ca parameters^d
	Na2O	CaO	B2O3	SiO2	P2O5	x NBO	x ^[4]B(MD)	x ^[4]B(NMR)	ε (%)	Ref.	Z̄ Na	Z̄ Ca	r̄(Na–O) (pm)	r̄(Ca–O) (pm)	x(Na–O^[1])	x(Ca–O^[1])
a Glass compositions modeled by atomistic MD simulations, with the sample notation described in Section 2 and details about the MD simulations provided in Table S1 (ESI). The NCPS and NCBPS glass compositions were taken from Mathew et al.^36,43 and Yu et al.,^30,44 respectively, while the NBS compositions were adapted from ref. 5, 8 and 21–23. b MD-derived NBO (O^[1]) fraction out of all O species, with the remaining constituting BO: xNBO + xBO = 1. c Fractional population x^[4]B of B^[4] coordinations obtained either by MD simulations or by ^11B NMR from the as-indicated literature source. The relative {B^[3], B^[4]} populations {x^[3]B, x^[4]B} obey x^[3]B + x^[4]B = 1. ε is the (signed) relative deviation between the modeled and experimental results: 100{x^[4]B(MD) − x^[4]B(NMR)}/x^[4]B(NMR). d Average coordination number (Z̄M), average M–O distance [r̄(M–O)], and fraction of O^[1] species [x(M–O^[1])] present in the ensemble of {MOp} polyhedra for M = {Na, Ca}. Note that x(M–O^[1]) + x(M–O^[2]) = 1. e The data listed in the upper and lower rows were obtained from slow and fast melt quenching, respectively. f The NCPS[2.54] base composition was used to design the NCBPS(f) glass series by progressive SiO2 → B2O3 replacements.^30,44 g Typical data uncertainties, estimated as the root-mean-square (rms) deviation from the average parameter value. The uncertainty of the NMR-derived x^[4]B data is that reported by Yu et al.,^44 which may be taken as representative also for the NB/NBS glasses (in the absence of stated data uncertainties in many of these reports). h Based on the simulation/experimental uncertainties.
NB0.11	10.0		90.0			0.000	0.118	0.116	1.7	4	7.30		258.9		0.00
NB0.20	17.0		83.0			0.000	0.214	0.225	−4.9	4	7.23		258.8		0.00
NB0.25	20.0		80.0			0.007	0.331	0.320	3.4	27	6.99		258.7		0.02
NB0.50	33.3		66.7			0.050	0.420	0.430	−2.3	4	6.91		257.6		0.09
NB0.67	40.0		60.0			0.139	0.420	0.438	−4.1	33	6.56		255.9		0.22
NB1.00	50.0		50.0			0.337	0.331				6.14		253.3		0.45
NB1.30	56.5		43.5			0.495	0.238				6.00		252.0		0.60
NBS2–0.25(33)	7.7		30.8	61.5		0.005	0.240	0.230	4.3	28	6.21		259.1		0.02
NBS2–0.50(33)	14.3		28.6	57.1		0.020	0.433	0.450	−3.8	4	6.42		258.5		0.06
NBS2–0.75(33)	20.0		26.7	53.3		0.050	0.563	0.581	−3.1	21	6.54		258.3		0.12
NBS2–2.50(33)	45.4		18.2	36.4		0.445	0.385	0.350	10.0	5	5.78		251.8		0.61
NBS0.5–0.50(67)	25.0		50.0	25.0		0.033	0.435	0.420	3.6	19	6.83		258.4		0.07
NBS5.10–1.31(16)	17.7		13.5	68.8		0.097	0.615	0.830	−25.9	8	6.04		256.9		0.24
NBS4–4.00(20)	44.5		11.0	44.5		0.484	0.367	0.380	−3.4	23	5.63		251.1		0.66
								0.330	11.2
NCBPS(0)^f	24.1	23.3		48.6	4.0	0.625					5.84	5.91	252.1	246.2	0.74	0.91
NCBPS(10)	24.1	23.3	4.9	43.7	4.0	0.584	0.403	0.388	3.9	30	5.95	6.03	252.4	247.1	0.70	0.87
NCBPS(20)	24.1	23.3	9.7	38.9	4.0	0.546	0.392	0.424	−7.5	30	6.09	6.16	253.0	248.4	0.65	0.84
NCBPS(30)	24.1	23.3	14.6	34.0	4.0	0.504	0.431	0.432	−0.2	30	6.24	6.19	253.7	248.7	0.60	0.80
NCBPS(50)	24.1	23.3	24.3	24.3	4.0	0.438	0.420	0.433	−3.0	30	6.37	6.43	254.2	251.0	0.52	0.73
NCBPS(80)	24.1	23.3	38.9	9.7	4.0	0.348	0.419	0.419	0.0	44	6.61	6.64	255.5	253.1	0.40	0.63
NCPS[2.11]	24.6	26.7		46.1	2.6	0.688					5.87	5.96	251.6	246.1	0.77	0.91
NCPS[2.30]	24.2	26.4		45.4	4.0	0.677					5.94	6.00	251.8	246.6	0.77	0.91
NCPS[2.54]^f	24.1	23.3		48.6	4.0	0.625					5.84	5.91	252.1	246.2	0.74	0.91
NCPS[2.74]	20.2	22.2		55.0	2.6	0.544					5.82	5.85	252.9	246.4	0.67	0.88
NCPS[2.93]	17.9	23.3		54.8	4.0	0.529					5.89	5.94	253.3	247.4	0.67	0.88
σ						0.001	0.007	0.01	2.0^h	44	0.02	0.02	0.1	0.2	0.01	0.01

---

The nomenclature of the Na-based borate [NBR] and borosilicate [NBSK–R(f)] glass members adopts the composition parametrization RNa₂O–B₂O₃–KSiO₂ of Yun–Dell–Bray–Xiao (YDBX)^4,5 with R = x(Na₂O)/x(B₂O₃) and K = x(SiO₂)/x(B₂O₃). The YDBX structural model for borate/borosilicate glasses is introduced briefly in Section 4, with details provided in ref. 4 and 5. The NBR series invokes an increasing Na₂O content (R) at K = 0. To sample a reasonably wide composition range with a small number of glasses, the NBS set encompasses 4 compositions with K = 2 and increasing R, together with 3 glasses of distinct {K, R} pairs.

The NCPS[N̄^Si_BO] phosphosilicate glasses are labelled by their silicate network connectivity (N̄^Si_BO), owing to its direct relevance for the bioactivity.^69,70,73 The NCBPS glass series was designed from the NCPS[2.54] base glass by replacing SiO₂ by B₂O₃ at unaltered Na₂O and CaO contents and a fixed sum x(SiO₂) + x(B₂O₃).^30,44 Consequently, the NCPS[2.54] and NCBPS(0) compositions are identical. To assist comparisons between the borosilicate and BPS compositions, the value f in each NCBPS(f) and NBK-R(f) label represents the percentage of B₂O₃ out of the total B₂O₃ and SiO₂ contents: f = 100x(B₂O₃)/[x(SiO₂) + x(B₂O₃)]. Note that besides comprising Ca, most NC(B)PS glass networks are markedly more fragmented than their NB(S) counterparts, as encoded by the NBO fraction out of the total O speciation in the structure: x_NBO = x_NBO/(x_NBO + x_BO), where x_NBO + x_BO = 1.

3 Computational methods and force-field developments

3.1 Molecular dynamics simulations

All glass models were generated by atomistic MD simulations using the DLPOLY4.08 program,^74,75 utilizing NVT ensembles in a cubic box with periodic boundary conditions, with the box size and the number of atoms selected to match the experimental glass composition and density. The borate and borosilicate glass models comprised ≈3500 atoms, whereas the P-doped NC(B)PS counterparts involved 6000–7000 atoms to ensure sufficient statistics of their minor P speciations (≤4 mol% P₂O₅). Here we only provide the common simulation procedures and parameters (which were also employed for the force-field developments); the details are given in Table S1 (ESI†).

Each melt-quench protocol started from a random atom configuration that was equilibrated for 100 ps at 3500 K, whereupon the temperature was reduced by 100 K every 10 ps (10 K ps⁻¹ nominal cooling rate) to 300 K. Then followed another equilibration for 200 ps, from which the last 150 ps were averaged to yield the structural data. This procedure was completed 2–4 times (with different initial atom configurations for each glass composition), from which the average value and uncertainty of each reported structural parameter were derived. The equations of motion were integrated in steps of 0.2 fs using the velocity Verlet integrator, while the temperature was controlled using a gentle stochastic thermostat with a 1.0 ps time constant and a Langevin friction constant of 1.0 ps⁻¹. These simulation conditions (e.g., system sizes and equilibration stages) provide well-converged and reliable structural parameters,^36,58 whereas the effects from the choice of cooling rate are examined in Sections 4.1 and 5.3.

All computations utilized a polarizable shell-model potential,^45,55,57 where each cation carries its full formal charge,⁴⁵ but the O²⁻ species are represented by core (O_C) and shell (O_S) portions with charges z_C = +0.8482e and z_S = −2.8482e, respectively (z_C + z_S = −2), and corresponding masses m_C = 15.7994u and m_S = 0.2000u, where e is the elementary charge and “u” is the atomic mass unit. The core–shell units are connected by a harmonic potential with the force constant 74.92 eV Å⁻².⁴⁵ The interaction energy of two atom/ion species α and β at a distance r_αβ was parametrized by the following extended Buckingham potential,

which accounted for all short-range O_S–O_S and cation–O_S pair interactions, and was evaluated out to r_αβ = 0.8 nm. Long-range 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⁻⁶. All {A_αβ, ρ_αβ, C_αβ, D_αβ} values are listed in Table 2. For brevity, we will henceforth drop the αβ subscripts whenever the identity of the atom-pair is obvious from the context.

Table 2 Interatomic potential parameters

α–β	Pair potential^a				Ref.
	A αβ (keV)	ρ αβ (Å)	C αβ (eV Å^6)	D αβ (eV Å^12)
a The pair potential involves a modified Buckingham term; see eqn (1). Full formal charges were employed for all cations, whereas the O^2− ions are represented by core (OC) and shell (OS) parts. b These parameters were first utilized in our recent report on BPS glass structures.^44 c The three-body potential is represented by a truncated harmonic function , which was used in our recent work^36,42,43 and involves a modified functional of that reported by Tilocca et al.^57 It was evaluated out to r = 250 pm.^75
Na–OS	56.465	0.1939	0	0	55
Ca–OS	2.152	0.3092	0.099	0	55
Si–OS	1.284	0.3205	10.662	0	45
B–OS	0.472	0.3350	0	5.4	44
P–OS	1.750	0.2900	0	0	44
OS–OS	22.764	0.1490	27.880	0	45
OC–OS	0	0	0	0.1	36

---

α–β–α	Three-body potential^c			Ref.
	k αβα (eV rad^−2)	ρ αβα (Å)	θ ^0αβα
OS–Si–OS	5.48	2.030	109.47°	36
OS–P–OS	2.74	2.030	109.47°	36

---

The average coordination number of E (Z̄_E) corresponds to the (mean) number of O atoms in its first coordination shell; it was determined by integrating the pair distribution function [PDF; g_E–O(r)] out to its first minimum,⁴⁷ which throughout all models corresponded to cutoff radii of 205 pm for {B, Si, P} and 312 pm for {Na, Ca}. Noteworthily, a minor variation (±20 pm) around each of these “most suitable” cutoff radii had insignificant bearing on the resulting Z̄_E value. Note that as opposed to the “most probable E–O distance” obtained from the maximum of the PDF, we report the average E–O distance (“bond length”) throughout, denoted by r̄(E–O).

3.2 P–O force field optimizations

The new P–O parameters resulted by empirically refining those of Tilocca⁵⁷ to reduce the discrepancies of the MD-generated fractional populations {xⁿ_P} of the {Qⁿ_P} species relative to experiments, while preserving the already accurate modeling of local structural features of the PO₄ tetrahedron.

As explained in Section 5, the MD-generated BO-bearing phosphate population x¹_P is consistently overestimated at the expense of x⁰_P. The P–NBO affinity may be increased by reducing the “softness” parameter ρ in eqn (1), which in practice was accomplished by locating the {A, ρ} pair that maximized x⁰_P in models of the NCPS[2.93] glass. The {A, ρ} values of Table 2 resulted from a mesh of ρ-values, with the accompanying parameter A obtained from the constraint of keeping the modeled P–O^[1] bond-length ≈154 pm for all {A, ρ} pairs. Further evaluations of NC(B)PS glasses verified a global improvement of the modeled {xⁿ_P} data, regardless of the precise glass composition and/or network polymerization degree (Section 5).

3.3 B–O force field optimizations

As opposed to the P–O parameters that resulted by adjusting those of ref. 57, the B–O parameters were developed from scratch. This is to our knowledge the first B–O interatomic potential option developed within the shell-model framework. Notably, the Si–O and P–O potentials invoke three-body terms (see Table 2) for preserving the tetrahedral geometry,^45,57,75 whereas the B–O force field should operate in their absence to permit conversions between triangular and tetrahedral geometries for an accurate modeling of the borate speciation. An increase (decrease) of the “softness” parameter ρ promotes a larger relative BO₄ (BO₃) population.

To avoid unrealistically short B–O (or more precisely B–O_S) distances from shell model MD simulations of the small and “soft” B atom, it was necessary to invoke a non-zero repulsive parameter (D) together with the {A, ρ, C} entities in eqn (1). After verifying that the parameter C had insignificant bearing on the simulation outcomes, we proceeded with C ≡ 0. Each value of the {A, ρ, D} triplet was deduced sequentially by the following three-step protocol:

(1) A {ρ, D} grid was generated over the ranges 0.2 ≤ ρ ≤ 0.4 and 0 ≤ D ≤ 20. For each {ρ, D} pair, the parameter A was obtained by energy minimizations on crystalline borates (B₂O₃; CaB₂O₄; NaBO₂; Na₃CaB₄O₁₀), whose structures feature different {x^[3]_B, x^[4]_B} contributions. These optimizations were performed with the GULP program⁷⁶ and the procedure outlined in the ESI,† except that all lattice parameters were fixed in the present calculations, while the atom positions were allowed to vary freely in the presence of the force fields.

(2) For each value of ρ, D was selected from the {A, ρ, D} triplet that minimized the energy in stage (1).

(3) The remaining parameter ρ was subsequently deduced from MD-generated models of the NCBPS(50) glass to best match its NMR-derived {x^[3]_B, x^[4]_B} fractions (Table 1).

Further assessments on B-bearing crystalline phases (see the ESI†) and glasses (Section 4) confirmed that the as-obtained B–O force field offers accurate borate speciations, thereby evidencing a very good transferability among different glass systems and compositions.

4 Borate speciations in multicomponent glasses

Here we evaluate the predictive power of the new B–O potential to quantitatively reproduce the fractional populations {x^[3]_B, x^[4]_B} of {BO₃, BO₄} groups in models of B-bearing glasses of increasing complexity from the Na₂O–B₂O₃, Na₂O–B₂O₃–SiO₂, and Na₂O–CaO–B₂O₃–SiO₂–P₂O₅ systems (see Table 1). Note that x^[4]_B is often denoted as “N₄” in the literature. Fig. 1 plots the x^[4]_B results obtained either by MD simulations or by ¹¹B MAS NMR experiments against the parameter R.^4,5 The NMR data were reproduced from various sources (see the caption of Fig. 1), with each value closest to the respective simulation outcome indicated by a red circle.

Fig. 1 Fractional populations (x^[4]_B) of BO₄ groups obtained either by MD simulations (squares) or by NMR experiments (circles) and plotted against R = x(Na₂O)/x(B₂O₃) for (a) Na₂O–B₂O₃ and (b) Na₂O–B₂O₃–SiO₂ glasses. The cyan solid lines represent x^[4]_B data predicted by the Yun–Dell–Bray–Xiao model, whereas the gray rectangles indicate its associated regions “I–IV”.^4,5 The NBS glasses involve distinct values of K = x(SiO₂)/x(B₂O₃), as identified in (b), whereas all NB glasses correspond to K = 0. The experimental data were selected from various literature sources; for borates: ref. 4, 27, 29, 33, 38, and 81; for borosilicates: ref. 4–6, 8, 20–23, 28, and 32. Several independent literature reports are available for some glass compositions: each solid circle marks the most representative NMR-derived value [“NMR(ref)”], which is compared with the MD-generated data in Table 1. The red bar associated with the NBS4–4(20) glass depicts the range of NMR-derived BO₄ fractions obtained from glasses prepared with distinct melt-cooling rates.²³ Here and elsewhere: data uncertainties within/around the symbol sizes are not displayed.

For increasing Na₂O content of the NB/NBS glass, Fig. 1 reproduces the well-known trend of an initial BO₃ → BO₄ conversion up to the R-value for which x^[4]_B is maximized, x^[4]_B(max) = R_max, which according to the YDBX model occurs at R_max = 0.5 and R_max = (K/8 + 1)/2 for the NB and NBS systems, respectively.^4,5 In Fig. 1, the YDBX predictions are depicted by lines in cyan color. Note that the relative BO₄ population grows concurrently with the SiO₂ content of the borosilicate glass, encoded by the parameter K = x(SiO₂)/x(B₂O₃), and that x^[4]_B(max) is consistently higher than the corresponding borate analog (K = 0).^4,5 When the Na₂O content (R) of an NB glass is increased further, BO₄ groups progressively convert into NBO-bearing BO₃ moieties (Fig. 1a), which marks the onset of NBO formation.^4,5 In contrast, x^[4]_B remains close to x^[4]_B(max) throughout “region III” of the NBS system (Fig. 1b): according to the YDBX model, only Si–NBO contacts exist in structures conforming to region III, while both SiO₄ and BO₃ groups (but not BO₄) accommodate NBO species in regime IV.^4,5 While the MD-generated x_NBO values of Table 1 accord overall well with the YDBX prediction (not shown), minor NBO populations occur already for R < R_max in both NB and NBS glasses, as also reported previously.^{10,19,20,37,51–53}

Fig. 1 and Table 1 evidence an excellent accordance between the modeled and NMR-derived borate speciations, with the respective x^[4]_B values typically agreeing within the experimental/modeled data uncertainties, and well within the spread of NMR estimates originating from distinct studies of a given glass composition; see Fig. 1. Across all 18 B-based glasses compared (also encompassing the NCBPS system) in Table 1, the modeled BO₄ fraction typically reproduces the experimental counterpart within 95%. The largest discrepancies of the simulated x^[4]_B data generally result for the Si-richest NBS/NCBPS glasses (f ≲ 20%), with the NBS5.10–1.31(16) glass model featuring the globally largest relative difference (26%). The borate speciations of the NCBPS glasses were evaluated further in ref. 44 for an expanded glass ensemble.

We conclude that despite non-negligible discrepancies of the {x^[3]_B, x^[4]_B} values relative to experiments observed for a few Si-rich NBS/NCBPS glasses (f ≲ 20%), our B–O potential provides overall excellent predictions of the {BO₃, BO₄} speciations in both borate and boro(phospho)silicate glasses over wide composition ranges: the predictive power typically matches that observed from system-dedicated force fields involving glass-composition dependent B–O potential parameters to enhance the performance.^51–54 The present force field also appears to provide at least as accurate predictions as that recently introduced by Pacaud et al. for NBS glasses,¹⁰ which also involves constant B–O parameters that consider polarization effects. Assessments of five Na₂O–B₂O₃–SiO₂ glasses with variable compositions¹⁰ (yet across a narrower span than the present NBS glasses) revealed a 0–26% spread of relative deviations between the simulated and NMR-derived BO₄ populations, with average and median deviations of 11% and 8.5%, respectively.¹⁰

4.1 Cooling-rate trends

The temperature-dependent borate speciation of the melt induces a strong relationship between the cooling-rate and the {x^[3]_B, x^[4]_B} fractions of the resulting glass structure,^{9,10,23,24,26} particularly considering the ∼10 orders of magnitude more rapid simulated quenching relative to the experimental glass production. To gauge the melt-cooling behavior for our present simulation conditions and B–O potential, we performed calculations with variable quench rates (q) for the NBS2–0.25(33), NBS2–0.75(33), and NCBPS(30) glasses, whose {x^[4]_B} results are shown in Fig. 2. The precise cooling-rate responses are clearly composition-dependent, but all glasses manifest a general trend of elevated BO₄ populations for decreasing q, as expected.

Fig. 2 Cooling-rate (q) dependence of the MD-generated relative BO₄ fraction for the as-indicated NBS and NCBPS glasses. The vertical line marks the rate q = 10 K ps⁻¹ employed in all other simulations and corresponding to the {x^[4]_B} data listed in Table 1. Each of the three solid horizontal lines marks the NMR-derived experimental value, with the corresponding grey area indicating the experimental uncertainty (±σ), whereas the dotted horizontal lines show the predictions from the YDBX model.^4,5

The NBS2–0.75 composition reveals the largest variation between the highest cooling rate (100 K ps⁻¹) and that of 10 K ps⁻¹ utilized in all other simulations herein. Yet, when q is decreased further, the x^[4]_B value converges to the experimental result within its ±σ uncertainty (Fig. 2). Moreover, it is gratifying that the BO₄ population of the NBS2–0.25 counterpart remains close to the experimental value throughout the entire evaluated range 1 ≤ q/(K ps⁻¹) ≤ 100. The NCBPS(30) composition reveals an intermediate behavior of a seemingly continuous decrease of the BO₄ population when the cooling rate varies from 2 K ps⁻¹ to 100 K ps⁻¹, while the B-richer NCBPS(50) composition manifests much lower variations (see Fig. S1, ESI†).

Notwithstanding the composition-dependent effects, we conclude that usage of q = 10 K ps⁻¹—for which the B–O potential parameters were also optimized—combines reasonable computation times with accurately modeled borate speciations.

5 Phosphate speciations in (boro)phosphosilicate glasses

5.1 Bioactive phosphosilicate glasses

Given the beneficial properties of phosphate incorporation into Na₂O–CaO–SiO₂ glasses for their in vitro apatite formation,^71–73,77 assessments of composition–structure–bioactivity trends by MD modeling require accurate predictions of the {Qⁿ_P} speciation over the network-connectivity span 2.0 ≲ N̄^Si_BO ≲ 3.0 and particularly for glasses across the range 2.1 ≤ N̄^Si_BO ≤ 2.6 with ≤6 mol% P₂O₅ that is relevant for all attainable bioactive glass compositions.⁷³ For this regime, the Q⁰_P fractional population (x⁰_P) decreases for increasing glass network polymerization, but is independent of its P content.³⁶ However, current force fields used in classical MD simulations cannot accurately reproduce the markedly stronger propensity for P–NBO contacts relative to Si–NBO: this leads to overestimated Q¹_P (and Q²_P) populations at the expense of the orthophosphate counterpart that is experimentally evidenced to dominate the phosphate speciations in all NCPS glasses with sufficiently large Na⁺/Ca²⁺ reservoirs for balancing all negative charges of the Q⁰_P tetrahedra.³⁶

For glasses with variable silicate network connectivities, Fig. 3a contrasts the orthophosphate fraction predicted from MD simulations with that estimated by ³¹P MAS NMR in ref. 36. We first consider the widely studied BG composition 24.6Na₂O–26.7CaO–46.1SiO₂–2.6P₂O₅ of Hench, “45S5”,⁶⁷ which features N̄^Si_BO = 2.11 and an NMR-derived phosphate speciation comprising 96 ± 1% of Q⁰_P groups, while Q¹_P accounts for the remaining, i.e., {x⁰_P, x¹_P} = {0.96, 0.04}.³⁶ MD simulations with the new P–O force field offer an excellent prediction (x⁰_P = 0.93 ± 0.015; circle in Fig. 3a), which translates into a 3% relative discrepancy, which is almost within ±σ of the experimental/simulation data uncertainties (and well within ±2σ). In contrast, a markedly larger deviation (11%) resulted when using the P–O parameters of ref. 57 in an otherwise identical computation (x⁰_P = 0.85 ± 0.015; solid triangle in Fig. 3a). These values may be contrasted with other shell-model MD simulations of the 45S5 glass utilizing the P–O force field of ref. 57 and a melt-cooling rate of 10 K ps⁻¹, yielding x⁰_P values of 0.65 (ref. 57), 0.73 (ref. 34), 0.80 (ref. 61), and 0.82 (ref. 62), while calculations employing rigid-ion potentials generally give markedly worse estimates of x⁰_P ≈ 0.5.^59,62 The relatively large spread in simulation outcomes when using the same P–O force-field parameters partially stems from system-size variations, where x⁰_P < 0.80 generally resulted from (too) small ensembles comprising <2000 atoms in the simulation box. However, other subtle factors may also contribute, as discussed below.

Fig. 3 Orthophosphate fractions obtained by MD simulations or NMR experiments^30,44 for Na₂O–CaO–SiO₂–P₂O₅ glasses with variable network condensation degrees; the data are plotted against the (a) silicate network connectivity (N̄^Si_BO), or the (b) NBO fraction out of all O species (x_NBO). The plot in (b) also includes the results from the NCBPS glasses (reproduced from ref. 44). All modeled data in (b), as well as those labeled “new P–O” in (a), resulted from the P–O force field herein. The data with solid/open triangles were obtained by employing the P–O potential parameters of ref. 57 with otherwise identical simulation conditions (solid triangles; “previous P–O”), or with the conditions of Mathew et al.³⁶ (open triangles). The experimental and simulation uncertainties are σ ≈ 0.01 and σ ≈ 0.025, respectively. Except for the (new) data on NCPS[2.30], the experimental results from the NCPS and NCBPS glasses were reproduced from Mathew et al.³⁶ and Yu et al.,^30,44 respectively.

For increasing silicate network connectivity, Fig. 3a reveals progressively larger deviations between the modeled and experimental {Qⁿ_P} speciations that implies substantial discrepancies for the most condensed NCPS glass networks with N̄^Si_BO ≥ 2.7. Yet, it is gratifying that significantly better predictions result with the new P–O potential throughout the N̄^Si_BO-range considered in Fig. 3a. For the NCPS[2.93] glass with an experimental orthophosphate fraction of 0.81, MD-derived counterparts of 0.54 and 0.67 were observed with the “previous”⁵⁷ and new P–O force fields, respectively, which correspond to relative deviations of 33% and 17% to the experiment and an error reduction by ≈50% for the new force field.

For sufficiently large atom (≳5000) ensembles and melt-cooling rates ≲10 K ps⁻¹, the converged simulation outcome would be expected to depend predominantly on the choice of force-field parameters. However, other factors may also affect the results, as illustrated by the two data-sets represented by solid/open triangles in Fig. 3a: both were obtained using the same P–O potential parameters (see ref. 57), where the {x⁰_P} results shown by solid triangles employed identical simulation conditions as those described in Section 3.1, while the obviously different data set labeled “Mathew2014” was obtained under identical conditions,³⁶ except for employing a Berendsen thermostat and an older version (3.10.0) of the DLPOLY program. Further evaluations (not shown) suggested that the distinct DLPOLY versions accounted for most of the large discrepancies between the two data sets, which are far outside of the statistical uncertainties, and almost as large as between the “new” and “previous”⁵⁷ P–O force fields in Fig. 3a.

5.2 Borophosphosilicate glasses

The significant formation of Qⁿ_P moieties with n > 0 in the glass models stems from difficulties of P–O force fields to reproduce the very strong P–NBO affinity in dense phosphosilicate networks (i.e., with low NBO contents) at practically accessible melt-cooling rates of ≳1 K ps⁻¹. This is illustrated by Fig. 3b, which plots the orthophosphate population against the NBO fraction of Na₂O–CaO–(B₂O₃)–SiO₂–P₂O₅ glasses. The equimolar B₂O₃-for-SiO₂ substitution in the NCBPS glass design (Section 2) is accompanied by a significant network condensation for increasing B₂O₃ content toward the limiting Na₂O–CaO–B₂O₃–P₂O₅ composition,⁴⁴ because the number of network-forming atoms doubles (1Si → 2B) while the corresponding number of O atoms only grows by 1.5 (2O → 3O).

For the most fragmented glass networks associated with x_NBO ≳ 0.55, Fig. 3b manifests modest deviations between the experimental and modeled x⁰_P values, as well as similar results between the NCPS and NCBPS glasses. However, as the NCBPS networks repolymerize (i.e., x_NBO is decreased), the modeled Q⁰_P populations diminish markedly, while the experimental counterparts alter marginally. The comparatively larger discrepancies between the NMR/MD-derived {x⁰_P} data observed for the NCBPS glasses relative to their NCPS counterparts for x_NBO ≈ 0.55 are attributed to the additional competition for NBO accommodation from a third network former (boron), which accentuates the underestimation of the P–NBO contacts; see ref. 44 and Section 6.1.

We conclude that relative to previous options, our new P–O force field greatly improves the modeled phosphate populations across the ranges 2.0 ≤ N̄^Si_BO ≤ 2.9 and x(P₂O₅) ≤ 0.06 that encompass the entire region of the Na₂O–CaO–SiO₂–P₂O₅ system relevant for rationalizing structure–composition–bioactivity trends. Yet, despite that the force field also offers substantial improvements for NCBPS glasses, non-negligible deviations between experiments and modeled {xⁿ_P} data remain that reflect limitations of (current) pair-potentials in classical MD simulations to accurately model the relative affinities among distinct glass-network formers to accommodate the BO/NBO species.

5.3 Cooling-rate trends

Analogous with the cooling-rate dependence of the borate speciation, the underestimated Q⁰_P populations observed in MD simulations stem partially from the orders-of-magnitude faster melt quenching necessitated in the calculations compared with the most rapid cooling attainable experimentally. Tilocca⁵⁸ examined the dependence of the modeled {Qⁿ_P} speciation of the 45S5 glass composition on the quench rate by using the P–O potential parameters of ref. 57, inferring that a limiting/converged value of x⁰_P ≈ 0.10 is expected for quench-rates ≲2 K ps⁻¹. Notably, this is more than twice that (x⁰_P = 0.04) of the hitherto sole experimental (non-zero) estimate³⁶ (note that ref. 58 involved comparisons with the experimental result x⁰_P = 0.08 reported from a Na-free 45S5 analog⁴⁰). The results herein demonstrate that MD-derived phosphate speciations closer to experiments are available by instead improving the P–O interatomic potential parameters while keeping the more attractive melt-cooling rate ≈10 K ps⁻¹ that allows for much faster computations.

We evaluated variable quench-rate MD simulations for the more difficult NCBPS scenario, where the improved P–O force field still yields markedly underestimated Q⁰_P populations. Fig. S1 (ESI†) reveals the strongest increase in the orthophosphate fraction when the quench-rate is reduced from 100 K ps⁻¹ to 16 K ps⁻¹. The growth of x⁰_P proceeds in the regime 2.5 ≤ q/(K ps⁻¹) ≤ 16, but to a lower extent, suggesting that modest improvements are expected for the prohibitively time-consuming simulations associated with q < 1 K ps⁻¹.

6 Local environments of the network formers

6.1 NBO distribution among network formers

Yu et al.⁴⁴ discussed the relative affinity of each network former {B^[3], B^[4], Si, P} to coordinate NBO in NCBPS glasses, mainly concluding that the propensity for NBO-accommodation decreases along the series P^[4] ≫ B^[3] > Si > B^[4]. Here we examine these NBO-partitioning trends among the network formers further by also considering glasses from the NB and NBS systems, none of which comprise P, whose very substantial NBO affinity (with an average of 3.0–3.8 NBO ions per tetrahedron) consumes ≈30% of the entire NBO reservoir in the NCBPS glasses even at the modest amount of 4 mol% P₂O₅.⁴⁴

Fig. 4 plots the average number of NBO ions accommodated by each F = {B^[3], B^[4], Si} species against x_NBO, in the NCBPS scenario employing the NBO fraction effectively available for B and Si (x^eff_NBO), i.e., after subtracting the P-associated NBO portion. The previously concluded relative NBO affinities among the various network formers⁴⁴ also hold in the less complex NB and NBS systems: notwithstanding that only three O positions are available at the BO₃ triangles, their NBO accommodation is higher than that of SiO₄ throughout all glass compositions, including regions “I–III” (see Section 4), which is at odds with the YDBX prediction.^4,5 Yet, similar observations were made in previous modeling studies of NBS glasses,^10,51–53 also concerning the (unexpected) presence of B^[4]–NBO contacts.^10,51 Indeed, Fig. 4 reveals an average number of 0.25–0.8 NBO ions at the BO₄ tetrahedra in the NCBPS networks, but also non-negligible B^[4]–NBO contacts in the NB(S) glasses with 0.3 < x_NBO < 0.5, in contrast with their assumed (near) absence in “conventional” borosilicate glass models^5,25,78 across the x_NBO range considered herein.

Fig. 4 MD-derived average number of NBO ions per network-forming species F = {B^[3], B^[4], Si} plotted against x_NBO of glasses from the (a) NB, (b) NBS2–R(30), and (c) NCBPS glass series. Open symbols in (b) are results for NBS4–4.00(20), whereas the dotted rectangle marks the plot-range employed for the inset graph. The NBO consumption by the PO₄ groups in the NCBPS glasses is accounted for in (c) by plotting against x^eff_NBO (see Section 6.1).

When compared at similar NBO contents, there is an overall increase in the number of NBO ions at each F species when progressing along the NB → NBS → NCBPS glass systems (even when disregarding the NBO consumption by P in the NCBPS models). Typically, an NBS/NCBPS glass network with x_Si/x_B ≈ 1 and x_NBO ≈ 0.5 comprises primarily BO₃ groups with 1 and 2 NBO ions, comparable amounts of BO₄ tetrahedra with 0 and 1 NBO ions, while Q³_Si groups constitute ≈50% of the silicate speciation, together with similar Q⁴_Si and Q²_Si populations.

6.2 F–O distances

Fig. 5(a, b) shows representative pair distribution functions associated with each network former F = {B^[3], B^[4], Si, P} in NB, NBS, NCPS, and NCBPS glasses with comparable network polymerization (x_NBO ≈ 0.5). Since both NBO (O^[1]) and BO (O^[2]) species are present in the glass networks, each PDF shown in the top panel of Fig. 5 has contributions from the F–O^[1] and F–O^[2] component plotted in the mid/bottom panels, with the average F–O distance, r̄(F–O), being the weighted average over the two r̄(F–O^[1]) and r̄(F–O^[2]) components, according to the BO/NBO distribution in the first coordination sphere of F.

Fig. 5 MD-derived pair distribution functions (PDFs) involving O and each network-former F = {B^[3], B^[4], Si, P} in the NB1.3, NBS2–2.5(33), NCPS[2.93], and NCBPS(30) glasses (see legends), which were selected to yield comparable NBO fractions x_NBO = {0.50, 0.45, 0.53, 0.50}. The PDFs involve the (a, b) total F–O distributions, as well as the (c–f) underlying F–O^[1] and F–O^[2] component distributions. Dotted vertical lines mark the various PDF maxima, specified by the distance in pm, while Table 3 lists the accompanying bond lengths.

Owing to the high field-strength of each network-forming cation, it tightly controls its first coordination shell, which remains essentially unperturbed by the neighboring F sites. Indeed, provided that the BO/NBO speciation is constant (x_NBO is fixed), each average distance r̄(F–O), r̄(F–O^[1]), and r̄(F–O^[2]) is nearly identical for a given F = {P, Si, B^[3], B^[4]} species both across and among the glass systems; see Table 3. This is expected from the very similar PDF curves in Fig. 5. The average F–O^[1] distance amounts to ≈154 pm for P, 161–162 pm for Si, while those of B^[3] and B^[4] are almost equal (≈134–135 pm). For each network former, the F–O^[1] distance distribution is narrower and its associated bond-length is slightly shorter relative to the F–O^[2] counterparts, with the typical r̄(F–O^[2]) − r̄(F–O^[1]) difference being ≈1 pm for P, ≈4 pm for Si and B^[3], but significantly larger for B^[4] (10–12 pm). These trends match the sequence of decreasing cation field-strength [CFS; ion charge divided by r̄(F–O^[q])²]: P > B^[3] ≳ Si > B^[4].

Table 3 Cation–oxygen bond lengths^a

Glass	x NBO	P–O			Si–O			B–O
		P–O	P–O^[1]	P–O^[2]	Si–O	Si–O^[1]	Si–O^[2]	B–O	B^[3]–O	B^[3]–O^[1]	B^[3]–O^[2]	B^[4]–O	B^[4]–O^[1]	B^[4]–O^[2]
a Average F–O^[1] and F–O^[2] distances in pm, as well as the F–O counterpart representative for all O species coordinated by each network former F = {P, Si, B^[3], B^[4]}. Each average distance was obtained by scanning over all FOp polyhedra in the glass model, only accounting for O species within 205 pm from F. b Typical data uncertainties, estimated as the rms deviation from the average parameter value.
NB0.11	0.000							136.1	134.9		134.9	142.6		142.6
NB0.20	0.000							137.1	135.1		135.1	142.6		142.6
NB0.25	0.007							138.4	135.4	133.9	135.4	142.9	134.1	142.9
NB0.50	0.050							139.5	135.7	134.0	135.8	143.4	134.2	143.4
NB0.67	0.139							139.9	136.1	134.2	136.4	143.9	134.5	144.0
NB1.00	0.337							139.9	136.7	134.4	137.7	144.6	134.6	145.2
NB1.30	0.495							139.6	137.0	134.6	138.7	145.7	134.7	147.2
NBS2–0.25(33)	0.005				163.3	159.6	163.3	136.8	134.5	133.6	134.5	142.3		142.3
NBS2–0.50(33)	0.020				163.7	160.1	163.8	138.6	134.9	133.8	134.9	142.3	134.5	142.3
NBS2–0.75(33)	0.050				164.1	160.4	164.1	139.9	135.4	134.0	135.5	142.5	134.6	142.6
NBS2–2.50(33)	0.445				164.8	161.6	165.9	140.5	136.9	134.7	138.9	144.9	135.0	146.8
NBS0.5–0.50(67)	0.033				164.3	160.3	164.3	139.2	135.4	133.9	135.5	142.9	134.2	142.9
NBS5.10–3.31(16)	0.097				163.9	160.3	164.1	140.1	135.5	133.9	135.8	142.3	134.9	142.5
NBS4–4.00(20)	0.484				164.8	161.7	166.1	140.2	137.0	134.7	139.8	144.4	135.0	146.2
NCBPS(0)	0.625	153.7	153.6	155.6	164.7	162.1	166.2
NCBPS(10)	0.584	153.7	153.6	155.2	164.7	161.9	166.0	140.6	136.7	134.7	139.1	145.0	135.1	147.4
NCBPS(20)	0.546	153.7	153.6	154.8	164.7	161.7	165.8	140.6	136.7	134.6	138.9	145.1	135.0	147.3
NCBPS(30)	0.504	153.7	153.6	154.6	164.7	161.7	165.6	140.9	136.6	134.6	138.3	145.1	134.9	147.2
NCBPS(50)	0.438	153.7	153.6	154.6	164.8	161.4	165.4	140.6	136.6	134.4	137.9	144.7	134.8	145.9
NCBPS(80)	0.348	153.6	153.5	154.2	164.9	161.2	165.1	140.4	136.6	134.2	137.1	144.6	134.6	145.3
NCPS[2.11]	0.688	153.7	153.7	155.7	165.0	162.7	167.0
NCPS[2.30]	0.677	153.7	153.7	155.7	164.9	162.4	167.0
NCPS[2.54]	0.625	153.7	153.6	155.6	164.7	162.1	166.2
NCPS[2.74]	0.544	153.7	153.6	155.8	164.6	161.8	165.9
NCPS[2.93]	0.529	153.7	153.5	155.8	164.4	161.6	165.5
σ	0.001	0.02	0.02	0.2	0.03	0.06	0.04	0.07	0.03	0.03	0.07	0.08	0.08	0.1

---

We next consider the composition dependence of the F–O/O^[1]/O^[2] bond lengths across the set of models from each glass system. Regardless of the glass network composition, the highest-CFS P⁵⁺ cation manifests an essentially constant distance of r̄(P–O^[1]) ≈ 153.6 pm and a very weakly varying P–O^[2] mean distance spanning 155.0 ± 0.8 pm. The P–O bond-length of 153.7 pm accords well with the diffraction-derived counterpart (155 ± 2 pm) for the Q⁰_P groups in the 45S5 BG.¹⁸ Somewhat larger (yet minor) variations are observed within the {r̄(F–O^[1])} and {r̄(F–O^[2])} distance-sets for Si, B^[3], and B^[4] (that manifest the same trends): the precise BO/NBO speciation of the glass accounts for most variations,^12,35,49,60 with both F–O^[1] and F–O^[2] distances increasing concurrently with x_NBO (see Table 3).

Table 3 reveals that the trend of increasing F–O^[1]/O^[2] bond lengths (for Si, B^[3], and B^[4]) offsets the expected shortening of the F–O counterpart for increasing x_NBO (and accompanying larger contributions from the shorter F–O^[1] distances): the bond-length merely either remains essentially constant or increases slightly, such as the r̄(Si–O) lengthening observed across the NCPS series. Notably, diffraction studies have also observed a Si–O distance-lengthening for increasing x_NBO.^12,79 Whereas experimental reference data are sparse even for simple binary silicate/borate glasses, we conclude that our modeled Si–O distances (163–165 pm) are typically 2–3 pm longer than diffraction-derived results from SiO₂, binary silicate glasses, and the 45S5 NCPS composition.^11–14,18 Moreover, it is gratifying that the MD-derived B–O and B^[3]–O bond-lengths of 139.5 pm and 135.7 pm match very well the experimental counterparts of 139.5 ± 1 pm and 137 ± 1 ppm for the NB0.50 glass composition,¹⁵ while a slightly larger deviation is observed between the modeled and experimental B^[4]–O distances of 143.4 pm and 147 ± 1 ppm, respectively. Noteworthily, the results of Table 3 reveal for each glass composition a 7–8 pm longer B^[4]–O average distance than B^[3]–O, in good agreement with previous findings by experiments and MD simulations.^10,15,53,54

For all non-negligible F–O^[2]–F′ linkages, we also evaluated the average distances between the F–F′ species present in neighboring polyhedra, which tend to decrease in the following order:

Si–Si(308–311) > Si–P(297–303) > B^[4]–P(290–297) > Si–B^[4](276–282) ≳ Si–B^[3](275–278) > B^[4]–B^[4](251–258) ≈ B^[3]–B^[4](251–258) > B^[3]–B^[3](243–248), (2)

with the distance-span (in pm) observed across all glasses given within parentheses. Besides depending on the intertetrahedral bond angles, these longer-range distances follow the gross trends anticipated from their r̄(F–O^[2]) and r̄(F′–O^[2]) components in Table 3.

6.3 Bond-angle distributions

As for the F–O distances, the O–F–O intrapolyhedral bond-angle distributions (BADs) plotted in Fig. S2 (ESI†) for the {BO₄, BO₃, SiO₄, PO₄} groups of the NB1.3, NBS2–2.5(33), NCPS[2.93], and NCBPS(30) glass networks confirm that each coordination geometry remains invariant despite the presence of multiple network formers in the glass: narrow distributions are observed for all polyhedra, albeit with slightly broader distributions for the BO₄ and BO₃ moieties. Hence, each Si/P/B^[3]/B^[4] network former manifests well-defined coordination polyhedra with essentially fixed average O–F–O intratetrahedral angles of 109.4° for PO₄, 109.3° for SiO₄, and 109.1° ± 0.15° for BO₄, while the planar BO₃ groups exhibit the mean angle 119.7° ± 0.1°. These values agree well with both experimental^11,16,17 and modeled^10,52,60 bond angles in (boro)silicate glasses. The local geometries are independent of the precise glass composition (within 0.3°), which is unsurprising for the Si–O and P–O force fields that exploit three-body terms for maintaining a strict tetrahedral geometry, but the new B–O counterpart apparently accomplishes this task even without such precautions (see Section 3.3).

We next consider the interpolyhedral bond-angle distributions of the F–O–F′ linkages involving {Si, B^[3], B^[4]} in the B/Si-bearing glasses. Table S3 (ESI†) compiles the average angles [small theta, Greek, macron](F–O–F′), and Fig. 6 plots the BAD functions for a selection of glasses with comparable network polymerization (x_NBO ≈ 0.5). The Si–O–Si BAD function is shifted to higher bond angles than all other {F, F′} pairs, as also reflected by the large value of [small theta, Greek, macron](Si–O–Si) ≈ 140° (Fig. 6e), whereas the other linkages exhibit similar mean angles [small theta, Greek, macron](F–O–F′) in the range 124°–136°; see Table S3 (ESI†) and Fig. 6.

Fig. 6 Interpolyhedral F–O–F′ bond-angle distributions observed in borate [NB1.3], borosilicate [NBS2–2.5(33)], and borophosphosilicate [NCBPS(30)] glass models for the as-indicated F–O–F′ linkages among the Si, B^[3], and B^[4] network formers. The numbers at the bottom left portion of each graph represent the average bond-angles shown with the same color coding for each glass as in the legend of (b). The peak ≈90° in (f) stems from a minor degree of BO₄–BO₄ edge-sharing [excluded in the averages provided in (f)].

The [small theta, Greek, macron](F–O–F′) values of Table S3 (ESI†) reflect the well-known trend of decreasing bond angles for increasing NBO content of the glass,^2,24,49,60 with the Si–O–Si/B^[3] angles showing the largest variations (within ≈10°) across a given glass series. Naturally, the fragmented NCPS networks that accompany the high NBO contents manifest relatively low [small theta, Greek, macron](Si–O–Si) angles compared to SiO₂, where the most recent studies confine [small theta, Greek, macron](Si–O–Si) ≈ 144°–147°.² Experimental data on multi-component glasses are very sparse, but an NMR/DFT study reported [small theta, Greek, macron](Si–O–Si) ≈ 133° for a Na₂O–SiO₂ glass with x_NBO = 0.58,³⁵ which may be contrasted with the values ≈139° observed from our NCPS models with comparable NBO contents (Table S3, ESI†). The bond-angle dependence on x_NBO partially accounts for the seemingly larger deviations observed among the glasses for the B^[4]–O–B^[4] pair in Fig. 6f, along with higher uncertainties of the BAD function for these comparatively rare interconnectivities (Section 6.1). Also noteworthy is the small peak ≈90° observed solely for the θ(B^[4]–O–B^[4]) distribution of each glass, which reflects a minor fraction of edge-shared BO₄ tetrahedra.

7 Structural roles of Na and Ca

7.1 Local Na and Ca environments

Table 1 lists the MD-derived average coordination numbers for Na (Z̄_Na) and Ca (Z̄_Ca) in each glass structure. The phosphosilicate glasses reveal very similar values of ≈5.9 for both Z̄_Na and Z̄_Ca, as discussed by Mathew et al.⁴³ All B-based glasses manifest higher mean Na/Ca coordination numbers, both growing for increasing B content or decreasing x_NBO. However, the {x_B, x_NBO} parameters are usually coupled by the glass design (see Table 1): for instance, the decrease in x_NBO is most likely mainly responsible for the reduced {Z̄_Na, Z̄_Ca} values when x(B₂O₃) grows along the NCBPS series. Variable NBO contents presumably also underlie similar Z̄_Na trends reported earlier for NBS and NCBPS glasses for increasing amount of B₂O₃.^50,54 At a fixed B content, Z̄_Na also increases concomitantly with the BO₄ population, which is particularly evident for the most polymerized borate/borosilicate networks with x_NBO ≲ 0.05 in Table 1.

Depending on the precise glass composition, the underlying {M^[p]} ensembles may comprise coordination numbers 4 ≤ p ≤ 9 (see Table S4, ESI†), but the distribution-width σ_M(p) ≈ 1 for both Na and Ca constrains only three polyhedral types to dominate in each structure: MO₅, MO₆, and MO₇, where MO₆ and MO₇ are most abundant when Z̄_M ≈ 6 and Z̄_M ≈ 7, respectively. Noteworthily, all cations coordinate both O^[1] and O^[2] species. Yet, there is a strong propensity for NBO-accommodation at the NaO_p and CaO_p polyhedra,^34,43,48,57 as evident from their associated x(Na–O^[1]) and x(Ca–O^[1]) data listed in Table 1 that represents the NBO fraction coordinated by the respective {NaO_p} and {CaO_p} ensemble. These Na/Ca–O^[1] fractions are consistently higher than x_NBO, notably so for the higher-CFS Ca²⁺ ion, whose coordination shells involve 63–91% of NBO contacts in the NC(B)PS glass models, whereas the corresponding range for the fraction of Na–O^[1] bonds is 35–69%. The strong M–O^[1] bonding preference has bearing on several other structural features discussed below.

We next consider the Na–O and Ca–O average distances plotted against x_NBO in Fig. 7, which also includes the r̄(M–O^[1]) and r̄(M–O^[2]) data for M = {Na, Ca}. Selected PDFs are depicted in Fig. S3 (ESI†). For both cations, the M–O^[1] bond lengths are consistently shorter than M–O^[2]. Yet, whereas the mean Ca–O^[1] distance is shorter than that of Na–O^[1] (by ≈5 pm) for a given glass composition, the situation is reversed for the O^[2] contacts, where Ca manifests (on the average) ≈12 pm longer distances than Na: throughout all NC(B)PS glasses, well-confined values of r̄(Na–O^[1]) ≈ 248 pm and r̄(Ca–O^[1]) ≈ 243 pm are observed, whereas the respective r̄(M–O^[2]) distance-spreads are larger, amounting to 260–266 pm for Na and 270–281 pm for Ca. However, notwithstanding that the various glass systems reveal similar spans of {r̄(Na–O^[2])}, the Na–O^[1] bond-length tends to increase concurrently with x_NBO in the NB and NBS systems, corresponding to the range 238–246 pm in Fig. 7b.

Fig. 7 Average M–O distances plotted against x_NBO for (a–c) Na and (d–f) Ca to each as-indicated O/O^[1]/O^[2] species and glass system; see the legend in (a). The solid and open symbols in (a–c) represent data for NBS glasses with K = 2 and K ≠ 2, respectively. The span of the vertical plot-range (48 pm) is around 50% of the full Na–O and Ca–O distance-spreads; see Fig. S3 (ESI†) for a selection of PDFs.

The net M–O bond lengths plotted in Fig. 7(a, d) evidence overall slightly shorter r̄(Ca–O) values (246–253 pm) compared to r̄(Na–O), the latter spanning 252–255 pm for the mixed-cation NC(B)PS glasses and 251–259 pm across all Na-bearing glasses. The average distances tend to shorten when x_NBO is increased, as expected from the concurrently growing number of contributing shorter M–O^[1] bonds (relative to M–O^[2]) and the preference for Na, and particularly Ca, to coordinate NBO rather than BO species (in contrast with the results for the F–O distances; see Section 6.2). These trends also accord with the well-known decrease in r̄(M–O) for a concomitant reduction in Z̄_M. Note that we report the average Na–O (and Ca–O) distance, as opposed to the “most probable” counterpart observed at the PDF maximum, which appears to be most frequently reported for the “bond length” in the literature (see for instance ref. 13, 15, 48, 54, 57, 60, and 61). Our modeled M–O distances agree well with both—meaning that the PDF maxima around 235–240 pm and 230–235 pm of the respective Na–O and Ca–O PDFs in Fig. S3 (ESI†) accord with those of ref. 18, 48, 54, 57, 60 and 61, while the respective average distances also agree well with literature data.^34,35

7.2 Modifier or compensator roles of Na⁺/Ca²⁺?

The commonly assumed roles of the M⁺/M²⁺ cations as either NBO-associated “network modifiers” or “charge compensators” of [BO₄]⁻ tetrahedra in B-based glasses^4–7,9,10 provide an intuitive and conceptually appealing qualitative structural model. Yet, experimental assessments of the alkali/alkaline-earth partitioning among the various anionic moieties are sparse and qualitative.^80,81 Inarguably, the M⁺/M²⁺ cations drive BO₃ → BO₄ conversions (Section 4.1) and indeed associate with [BO₄]⁻ tetrahedra, as proposed by the YDBX model^4,5 and also by earlier descriptions of borate/borosilicate glasses.^3,78 However, the categorical “modifier”/“compensator” view is oversimplified in (at least) two aspects: (i) despite that SiO₄ and BO₃ groups devoid of NBO ions are formally uncharged, their BO atoms are well-known to participate in the M⁺/M²⁺ coordination shells, as discussed above. (ii) A minor but x_NBO-dependent fraction of the BO₄ ensemble involves B^[4]–O^[1] contacts⁴⁴ (see Fig. 4), whose associated M⁺/M²⁺ cations feature a dual compensator/modifier role.

For every glass model, we assessed the partitioning of each M = {Na, Ca} species among the following four structural fragments, where M may act as (A) a “modifier”, encompassing all M–O^[1] contacts but those involving BO₄ groups with at least one B^[4]–O^[1] bond, or as (B) a “compensator” of BO₄ moieties devoid of NBO; M may also associate with (C) NBO-bearing borate tetrahedra, or (D) BO species associated with the BO₃ and SiO₄ moieties (but excluding BO₄). The relative contributions, {x_M(O^[1]), x_M(BO₄), x_M(BO₄–O^[1]), x_M(O^[2])}, associated with the corresponding structural scenarios A–D were determined by scanning over the {MO_p} ensemble, and assigning each of the p O sites at the MO_p moiety to its relevant category A–D (weighted by 1/p, such that summation over all M–O contacts in the structure only counts each M site once).

First considering Na, Fig. 8a–c illustrates the dependence of the x_Na(O^[1]), x_Na(O^[2]), x_Na(BO₄), and x_Na(BO₄–O^[1]) fractions on the NBO content in each B-bearing glass. These fractions are foremost dictated by the network polymerization degree and secondly by the relative Si/B contents, which for the NB/NBS glasses are parametrized by R and K, respectively (Table 1). All glass systems manifest the following gross trends:

Fig. 8 Dependence of the fractions {x_M(O^[1]), x_M(O^[2]), x_M(BO₄), x_M(BO₄–O^[1])} on the NBO content in (a) NB, (b) NBS, and (c, d) NCBPS glass models, with the results for (a–c) M = Na and (d) M = Ca. The solid and open symbols in (b) are data for the NBS2–R(33) series and the NBS4–4.00(20) composition, respectively.

(i) The Na ensemble mainly partitions between scenarios A and D, i.e., Na coordinates O^[1] and O^[2] species, respectively, implying that the x_Na(O^[1]) and x_Na(O^[2]) fractions typically dominate over all others in glasses with high (x_NBO ≳ 0.35) and low (x_NBO ≲ 0.35) NBO contents, respectively.

(ii) The “charge-compensator” role of Na is comparatively minor throughout: x_Na(BO₄) naturally correlates with the number of BO₄ moieties in the glass (i.e., with the product x_Bx^[4]_B), meaning that high values of x_Na(BO₄) ≳ 0.2 are only observed for B-rich and highly polymerized glass networks adhering to the YDBX regimes I–III, for which x^[4]_B is large;^4,5 see Section 4. However, even when x_Na(BO₄) is nearly maximized, the number of Na–O^[2] contacts involving BO₃/SiO₄ groups generally outnumbers those of Na–BO₄. Similar fractions of x_Na(BO₄) ≈ x_Na(O^[2]) are only observed for NB glasses with 0.5 ≲ R ≲ 0.7 (Fig. 8a), whereas x_Na(BO₄) is consistently lower than x_Na(O^[2]) for any borosilicate glass: the dependence of x_Na(BO₄) on the relative B/Si content (i.e., the parameter K) may be gauged by contrasting the higher x_Na(BO₄) values observed from the borate glasses (K = 0) in Fig. 8a with the (lower) counterparts from the NBS glasses with K = 2 in Fig. 8b. Moreover, the fraction of Na–O^[1] contacts grows rapidly when x_NBO is increased, such that Na⁺ ions coordinate both O^[1] (scenario A) and O^[2] (scenario D) species to a higher extent than associating with BO₄ groups. This feature is particularly evident for the NCBPS glasses with x_NBO ≳ 0.5 (Fig. 8c), but also applies for all NB(S) networks conforming to the YDBX region IV.

(iii) The minor but non-negligible number of B^[4]–O^[1] contacts (see Section 6.1 and ref. 44) implies that few Na⁺ species assume the dual modifier/compensator role (i.e., scenario C) throughout all glass structures, where typically x_Na(BO₄–O^[1]) < 0.05. Yet, the number of B^[4]–O^[1] contacts grows concurrently with each of x_NBO and x_Bx^[4]_B: naturally, the largest x_Na(BO₄–O^[1]) values are encountered in NBO-rich networks that simultaneously exhibit high BO₄ populations, such as in all NCBPS structures, as well as in the NB1.30 and NBS2–2.50(33) glasses [x_Na(BO₄–O^[1]) ≈ 0.07]; see Fig. 8a–c.

All qualitative trends in the Na⁺ partitioning among the four scenarios A–D also apply to the divalent Ca²⁺ cation in the mixed-modifier NCBPS glasses (Fig. 8d), with the primary distinction of further emphasized Ca–O^[1] contacts (see Section 7.1) at the expense of Ca–O^[2] and Ca–BO₄. Consequently, for a fixed glass composition, the x_Ca(O^[1])–x_Ca(O^[2]) difference observed in Fig. 8d is markedly larger than the x_Na(O^[1])–x_Na(O^[2]) counterpart of Fig. 8c. Moreover, in all highly fragmented glass networks (x_NBO ≳ 0.45), Fig. 4c reveals that ≈50% of all BO₄ tetrahedra accommodate at least one NBO ion, which rationalizes why both Na/Ca populations of case Cexceed that of B, i.e., why x_M(BO₄–O^[1]) > x_M(BO₄); see Fig. 8(c and d).

We wish to underscore two aspects: first, our analysis based on the total number of M–O contacts of the {MO_p} ensemble highlights the oversimplification of the “compensator/modifier” roles, where particularly the ignored M–O^[2] contacts are substantial. Second, we also considered NBO-rich borate/borosilicate glass networks, which have received much less attention in the literature. However, once excluding the M–O^[2] contacts from the statistics and focussing on highly polymerized NB(S) glasses, it is evident from Fig. 8(a, b) that the “compensator” portion of the {Na⁺} ensemble dominates its “modifier” counterpart in NB structures with x_NBO ≲ 0.2 and in NBS2–R networks with x_NBO ≲ 0.15, as confirmed experimentally for NB glasses.⁸¹ Notably, the relative “compensator/modifier” contributions then accord well with similar inferences from NBS glass models by Pacaud et al.,¹⁰ which were obtained by a different approach.

7.3 Partitioning of Na⁺ and Ca²⁺ among the network formers

We next move the spotlight onto the partitioning of the {Na⁺} and {Ca²⁺} ensembles around each network former, F = {B^[3], B^[4], Si, Q⁰_P, Q^n>0_P}, where “Q^n>0_P” comprises all BO-bearing phosphate groups (mainly involving Q¹_P). These fractional populations are denoted by x_Na(F) and x_Ca(F), respectively. This Na/Ca partitioning is complementary to that examined in Section 7.2 in that all M–O^[1] contacts, which were previously split into x_M(O^[1]) and x_M(BO₄–O^[1]), are considered for all NBO-associated FO_p polyhedra, which involve primarily SiO₄, BO₃, and (to a lesser extent) BO₄. Also, no distinction is made between BO₄ groups with/without NBO ions, whose relative M contacts are grouped together into x_M(B^[4]).

Fig. 9a–c plots each fraction x_Na(F) against x_NBO in the NB, NBS, and NCBPS glass networks. In the binary Na₂O–B₂O₃ system (Fig. 9a), Na is shared exclusively between B^[3] and B^[4] coordinations, where x_Na(B^[4]) initially grows at the expense of x_Na(B^[3]) for increasing Na₂O content, as anticipated (Section 4). The Na–B^[3] (Na–B^[4]) associations are lowest (highest) when the BO₄ population is maximized (Fig. 9a), but x_Na(B^[3]) remains larger than x_Na(B^[4]) throughout all glass compositions, as expected from the results of Section 7.2. Fig. 9b reveals the same qualitative trends of the Na–B^[3]/B^[4] contacts in the NBS glasses, but now a significant portion of the Na⁺ cations also associate with Si, with the partitioning among each silicate (SiO₄) and borate (BO₃ and BO₄) ensemble roughly given by the x_Si/x_B ratio, as evident by contrasting the results from glasses with distinct values of K. A nearly proportional partitioning of Na among Si and B also holds for the NCBPS networks (Fig. 9c), but a non-negligible fraction of the {Na⁺} reservoir now resides around the phosphate groups (despite their modest content). Fig. 9d manifests the same gross trends of the Ca partitioning as that of Na, with the main distinction that the Ca–BO₃ associations are comparatively more abundant at the expense of Ca–SiO₄ and Ca–BO₄. The relatively larger number of Na–BO₄ contacts stems from the sole B^[4]–BO bonding at most BO₄ tetrahedra, combined with the better charge-matching of Na⁺–[BO₄]⁻ than Ca²⁺–[BO₄]⁻.

Fig. 9 Partitioning of the (a–c) Na⁺ and (d) Ca²⁺ ensembles around the network formers F = {B^[3], B^[4], Si, P} in the (a) NB, (b) NBS, and (c, d) NCBPS glass models. The labels “P(n = 0)” and “P(n > 0)” in (c) and (d) are the results for orthophosphate (Q⁰_P) and all BO-bearing (Q^n>0_P) groups, respectively. The solid symbols in (b) represent data for the NBS2–R(33) series, whereas the open symbols are those for the as-indicated K values.

The trends in the various {x_M(F)} fractions may be rationalized from the relative amount of each {B^[3], B^[4], Si, P} network species coupled with its propensity to associate with Na⁺/Ca²⁺, which depends primarily on the net negative charge of each FO_p moiety. The relative preference for each Na/Ca–F contact was assessed by calculating P_M(F) = x_M(F)/x_F, which is larger than unity for a preferential association of M with F: the higher the P_M(F)-value, the stronger the preference. While the set of precise {P_M(F)} values depends on the glass composition, notably so on x_NBO, the following general trends apply (data not shown): (i) P_M(Q⁰_P) ≫ P_M(Q^n>0_P) > P_M(F) for F = {B^[3], B^[4], Si}. If disregarding the strong preference for Na–PO₄ (and particularly) Ca–PO₄ associations, the following observations are noteworthy in the NB(S) systems: (ii) There is a prominent preference for Na–B^[4] contacts in highly polymerized networks (x_NBO ≲ 0.05), as reflected in P_Na(B^[4]) values in the range 1.2–1.8, while P_Na(B^[3]) ≈ 1 and 0.86 ≲ P_Na(Si) ≲ 1, meaning that there are fewer Na–Si contacts than predicted by a statistical distribution of Na around the {Qⁿ_Si} moieties. (iii) Yet, the P_Na(B^[4]) values diminish rapidly for increasing NBO content of the network, with P_Na(B^[4]) remaining lower than both P_Na(B^[3]) and P_Na(Si) whenever x_NBO > 0.1. Altogether, these features rationalize why the “BO₄ compensator” portion of the {Na⁺} reservoir is only significant in highly polymerized NB(S) glass networks (see Section 7.2).

7.4 Relative propensities for Na/Ca–F contacts in BPS glasses

We finally consider the relative propensities for Na⁺ and Ca²⁺ to associate with each network-forming species in the NCBPS glasses. This aspect was discussed in detail for (B-free) NCPS glasses by Mathew et al.,⁴³ to which we refer for details. The key property behind all trends—also underlying those discussed in Sections 7.2 and 7.3—is the attraction strength between the positive Na⁺/Ca²⁺ cations and the negatively charged {FO_p} network species, with more (less) charged FO_p moieties exhibiting stronger (weaker) contacts with both Na⁺ and Ca²⁺, but where the preference for Ca²⁺–FO_p associations grows with the number of NBO sites at FO_p, as follows:

(i) FO_p moieties devoid of NBO ions prefer coordination of Na⁺ rather than Ca²⁺, meaning that a larger number of Na⁺ cations are present in the second coordination shell of F than that suggested by a statistical distribution and the x_Na/x_Ca ratio. (ii) FO_p polyhedra with one NBO ion exhibit an essentially statistical Na⁺/Ca²⁺ distribution according to the x_Na/x_Ca ratio. (iii) All network groups featuring at least two NBO species strongly prefer to associate with Ca²⁺.

We remind that these trends are qualitative and the details depend on both x_NBO and the precise speciation of all FO_p polyhedra in the glass network. Two exceptions are noteworthy: first, BO₄ tetrahedra with one NBO species manifest only a very weak preference for Ca coordination, despite that their net negative charge is equivalent to a BO₃/SiO₄ group coordinating two NBO ions. Second, the Q¹_P groups in the NCBPS glasses prefer coordinating Ca over Na, in contrast with the pronounced Na–Q¹_P preference observed in the NCPS counterparts, as discussed in ref. 43.

8 Concluding remarks

We have explored the short-range structural features of multi-component glasses with increasing complexity, ranging from the B-based Na₂O–B₂O₃ and Na₂O–B₂O₃–SiO₂ systems that involve Na as the sole modifier, to the more complex Na/Ca-bearing four-component Na₂O–CaO–SiO₂–P₂O₅ and five-component Na₂O–CaO–B₂O₃–SiO₂–P₂O₅ systems, where the BPS glasses involve a large number of co-existing network-associated {BO₃, BO₄, Qⁿ_Si, Qⁿ_P} polyhedra. This is to our knowledge the first MD simulation report providing comparisons among four distinct glass systems (comprising up to three network formers) and involving samples together spanning a large range of compositions within each system.

Despite the complexities, we demonstrated that accurate glass models are available by employing established shell-model Si–O/Na–O/Ca–O potentials,^45,55,57 but new B–O and improved P–O force fields, validated against experimental data with overall gratifying results: (i) the modeled/experimental {x^[3]_B, x^[4]_B} fractions accorded very well, typically with <5% relative deviations, and with <10% discrepancies for all but one of 18 evaluated B-bearing glasses. (ii) The new P–O potential offers significantly improved predictions of the orthophosphate population in bioactive Na₂O–CaO–SiO₂–P₂O₅ glasses across the entire (relevant) network polymerization range 2.0 ≤ N̄^Si_BO ≤ 2.9, with relative deviations of 3–17% to experiments. Yet, the discrepancies grow rapidly for more condensed glass networks, such as for B-rich NCBPS glasses. (iii) Na/Ca–O and Si/B/P–O average distances, as well as intra/inter-polyhedral bond angles, agreed well with literature data. Last but not least, further successful validations of the B–O and P–O force fields against experimental constraints on the medium-range structures of NCBPS glasses are provided in ref. 44.

Owing to the high field-strength of each F = {B^[3], B^[4], Si, P} network-forming cation, its local coordination environments remain essentially unperturbed, regardless of the presence of other network formers in the glass structure. Indeed, all of the minor alterations of F–O/O^[1]/O^[2] bond-lengths and F–O–F′ interpolyhedral bond angles observed among the glasses are readily attributed to variations of the precise BO/NBO speciation rather than the particular glass system considered.

Concerning the local {NaO_p} and {CaO_p} environments, 5, 6, and 7 coordinations dominate each distribution, which shifts towards higher coordination numbers for decreasing NBO content of the glass. Across the entire ensemble of 25 modeled glass compositions, Na revealed average coordination numbers between 5.8 and 7.2, whereas the span for Ca was 5.8–6.6 in the Ca-bearing NC(B)PS glasses. While NBO ions dominate both coordination shells of Na⁺ and Ca²⁺, the preference for Ca²⁺–NBO contacts is stronger. These features have bearing on the partitioning of each {NaO_p} and {CaO_p} ensemble among the various silicate, borate, and phosphate groups. The relative Na–F and Ca–F associations depend primarily on the net negative charge of the FO_p polyhedron, with those devoid of NBO species manifesting a preference for Na⁺ in the mixed-modifier NCBPS glasses, while the modifiers are statistically distributed around FO_p groups with one F–NBO bond, and those featuring more than one F–NBO contact prefer to associate with Ca²⁺. Moreover, the pronounced Na/Ca–NBO preference implies that the portion of each {Na⁺} and {Ca²⁺} reservoir that charge-balances the [BO₄]⁻ tetrahedra is only significant in B-rich borate/boro(phospho)silicate glass networks with low NBO contents, while the cations predominantly act as “modifiers” in all glasses with non-negligible NBO contents. Yet, even when the Na–BO₄ associations are maximal in a boro(phospho)silicate glass, they are nevertheless outnumbered by the contacts between Na and the BO sites of SiO₄/BO₃ groups. Moreover, out of all Na⁺/Ca²⁺–BO₄ associations in NBO-rich glass networks, a substantial fraction represents NBO-bearing borate tetrahedra.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This work was funded by the Swedish Research Council (VR-NT 2014-4667) and the Carl Trygger Foundation (CTS 16:123), with computer resources provided by the Swedish National Infrastructure for Computing (SNIC at NSC; project 2016-1-493).

References

1. G. N. Greaves and S. Sen, Adv. Phys., 2007, 56, 1–166.
2. M. Edén, Annu. Rep. Prog. Chem., Sect. C: Phys. Chem., 2012, 108, 177–221.
3. A. C. Wright, Phys. Chem. Glasses: Eur. J. Glass Sci. Technol., Part B, 2010, 51, 1–39.
4. Y. H. Yun and P. J. Bray, J. Non-Cryst. Solids, 1978, 27, 363–380.
5. W. J. Dell, P. J. Bray and S. Z. Xiao, J. Non-Cryst. Solids, 1983, 58, 1–16.
6. A. Grandjean, M. Malki, V. Montouillout, F. Debruycker and D. Massiot, J. Non-Cryst. Solids, 2008, 354, 1664–1670.
7. D. Manara, A. Grandjean and D. R. Neuville, Am. Mineral., 2009, 94, 777–784.
8. X. Wu, R. E. Youngman and R. Dieckmann, J. Non-Cryst. Solids, 2013, 378, 168–176.
9. F. Michel, L. Cormier, P. Lombard, B. Beuneu, L. Galoisy and G. Calas, J. Non-Cryst. Solids, 2013, 379, 169–176.
10. F. Pacaud, J.-M. Delaye, T. Charpentier, L. Cormier and M. Salanne, J. Chem. Phys., 2017, 147, 161711.
11. R. L. Mozzi and B. E. Warren, J. Appl. Crystallogr., 1969, 2, 164–172.
12. M. Misawa, D. L. Price and K. Suzuki, J. Non-Cryst. Solids, 1980, 37, 85–97.
13. G. N. Greaves, A. Fontaine, P. Lagarde, D. Raoux and S. J. Gurman, Nature, 1981, 293, 611–616.
14. D. I. Grimley, A. C. Wright and R. N. Sinclair, J. Non-Cryst. Solids, 1990, 119, 49–64.
15. J. Swenson, L. Börjesson and W. S. Howells, Phys. Rev. B: Condens. Matter Mater. Phys., 1995, 52, 9310–9319.
16. R. L. Mozzi and B. E. Warren, J. Appl. Crystallogr., 1970, 3, 251–257.
17. H. F. Poulsen, J. Neuefiend, H.-B. Neumann, J. R. Schneider and M. D. Zeidler, J. Non-Cryst. Solids, 1995, 188, 63–74.
18. R. A. Martin, H. L. Twyman, G. J. Rees, E. R. Barney, R. M. Moss, J. M. Smith, R. G. Hill, G. Chibin, T. Charpentier, M. E. Smith, J. V. Hanna and R. J. Newport, J. Mater. Chem., 2012, 41, 2212–2223.
19. G. El-Damrawi and W. Müller-Warmuth, J. Non-Cryst. Solids, 1992, 146, 137–144.
20. R. Martens and W. Müller-Warmuth, J. Non-Cryst. Solids, 2000, 265, 167–175.
21. L.-S. Du and J. F. Stebbins, J. Non-Cryst. Solids, 2003, 315, 239–255.
22. L.-S. Du and J. F. Stebbins, J. Phys. Chem. B, 2003, 107, 10063–10076.
23. S. Sen, Z. Xu and J. F. Stebbins, J. Non-Cryst. Solids, 1998, 226, 29–40.
24. F. Angeli, O. Villain, S. Schuller, T. Charpentier, D. de Ligny, L. Bressel and L. Wondraczek, Phys. Rev. B: Condens. Matter Mater. Phys., 2012, 85, 054110.
25. M. M. Smedskjaer, J. C. Mauro, R. E. Youngman, C. L. Hogue, M. Potuzak and Y. Yue, J. Phys. Chem. B, 2011, 115, 12930–12946.
26. J. Wu and J. F. Stebbins, J. Am. Ceram. Soc., 2014, 97, 2794–2801.
27. J. D. Epping, H. Eckert, A. W. Imre and H. Mehrer, J. Non-Cryst. Solids, 2005, 351, 3521–3529.
28. S. Wegner, L. van Wüllen and G. Tricot, Solid State Sci., 2010, 12, 428–439.
29. Y. Kim and K. Morita, J. Non-Cryst. Solids, 2017, 471, 187–194.
30. Y. Yu and M. Edén, RSC Adv., 2016, 6, 101288–101303.
31. G. Tricot, Phys. Chem. Chem. Phys., 2016, 18, 26764–26770.
32. T. Nanba, M. Nishimura and Y. Miura, Geochim. Cosmochim. Acta, 2004, 68, 5103–5111.
33. D. Zielniok, C. Cramer and H. Eckert, Chem. Mater., 2007, 19, 3162–3170.
34. A. Pedone, T. Charpentier, G. Malavasi and M. M. Menziani, Chem. Mater., 2010, 22, 5644–5652.
35. F. Angeli, O. Villain, S. Schuller, S. Ispas and T. Charpentier, Geochim. Cosmochim. Acta, 2011, 75, 2453–2469.
36. R. Mathew, B. Stevensson, A. Tilocca and M. Edén, J. Phys. Chem. B, 2014, 118, 833–844.
37. P. Zhao, S. Kroeker and J. F. Stebbins, J. Non-Cryst. Solids, 2000, 276, 122–131.
38. J. F. Stebbins, P. Zhao and S. Kroeker, Solid State Nucl. Magn. Reson., 2000, 16, 9–19.
39. E. Leonova, I. Izquierdo-Barba, D. Arcos, A. López-Noriega, N. Hedin, M. Vallet-Regí and M. Edén, J. Phys. Chem. C, 2008, 112, 5552–5562.
40. F. Fayon, C. Duée, T. Poumeyrol, M. Allix and D. Massiot, J. Phys. Chem. C, 2013, 117, 2283–2288.
41. R. Mathew, C. Turdean-Ionescu, B. Stevensson, I. Izquierdo-Barba, A. García, D. Arcos, M. Vallet-Regí and M. Edén, Chem. Mater., 2013, 25, 1877–1885.
42. B. Stevensson, R. Mathew and M. Edén, J. Phys. Chem. B, 2014, 118, 8863–8876.
43. R. Mathew, B. Stevensson and M. Edén, J. Phys. Chem. B, 2015, 119, 5701–5715.
44. Y. Yu, B. Stevensson and M. Edén, J. Phys. Chem. B, 2017, 121, 9737–9752.
45. M. J. Sanders, M. Leslie and C. R. A. Catlow, J. Chem. Soc., Chem. Commun., 1984, 1271–1273.
46. Q. Xu, K. Kawamura and T. Yokokawa, J. Non-Cryst. Solids, 1988, 104, 261–272.
47. X. Yuan and A. N. Cormack, J. Non-Cryst. Solids, 2001, 283, 69–87.
48. A. N. Cormack and J. Du, J. Non-Cryst. Solids, 2001, 293-295, 283–289.
49. J. Du and A. N. Cormack, J. Non-Cryst. Solids, 2004, 349, 66–79.
50. F. Gou, G. N. Greaves, W. Smith and R. Winter, J. Non-Cryst. Solids, 2001, 293–295, 539–546.
51. M. Barlet, A. Kerrache, J.-M. Delaye and C. L. Rountree, J. Non-Cryst. Solids, 2013, 382, 32–44.
52. L.-H. Kieu, J.-M. Delaye, L. Cormier and C. Stolz, J. Non-Cryst. Solids, 2011, 357, 3313–3321.
53. H. Inoue, A. Masuno and Y. Watanabe, J. Phys. Chem. B, 2012, 116, 12325–12331.
54. X. Lu, L. Deng, P.-H. Kuo, M. Ren, I. Buterbaugh and J. Du, J. Mater. Sci., 2017, 52, 8793–8811.
55. A. Tilocca, N. H. de Leeuw and A. N. Cormack, Phys. Rev. B: Condens. Matter Mater. Phys., 2006, 73, 104209.
56. A. Tilocca and A. N. Cormack, J. Phys. Chem. B, 2007, 111, 14256–14264.
57. A. Tilocca, A. N. Cormack and N. H. de Leeuw, Chem. Mater., 2007, 19, 95–103.
58. A. Tilocca, J. Chem. Phys., 2013, 139, 114501.
59. Y. Xiang and J. Du, Chem. Mater., 2011, 23, 2703–2717.
60. A. Pedone, G. Malavasi, M. C. Menziani, A. N. Cormack and U. Segre, J. Phys. Chem. B, 2006, 110, 11780–11795.
61. A. Pedone, G. Malavasi and M. Menziani, J. Phys. Chem. C, 2009, 113, 15723–15730.
62. A. Tilocca, J. Chem. Phys., 2008, 129, 084504.
63. J. K. Maranas, Y. Chen, D. K. Stillinger and F. H. Stillinger, J. Chem. Phys., 2001, 115, 6578–6589.
64. O. L. G. Alderman, G. Ferlat, A. Baroni, M. Salanne, M. Micoulat, C. J. Benmore, A. Lin, A. Tamalonis and J. K. R. Weber, J. Phys.: Condens. Matter, 2015, 27, 455104.
65. G. Ferlat, The B₂O₃ Case, in Molecular Dynamics Simulations of Disordered Materials, Springer International Publishing, 2015.
66. A. Zeidler, K. Wezka, D. A. J. Whittaker, P. S. Salmon, A. Baroni, S. Klotz, H. E. Fischer, M. C. Wilding, C. L. Bull, M. G. Tucker, M. Salanne, G. Ferlat and M. Micoulat, Phys. Rev. B: Condens. Matter Mater. Phys., 2014, 90, 024206.
67. L. L. Hench, J. Am. Ceram. Soc., 1991, 74, 1487–1510.
68. J. R. Jones, Acta Biomater., 2013, 9, 4457–4486.
69. Z. Strnad, Biomaterials, 1992, 13, 317–321.
70. R. Hill, J. Mater. Sci. Lett., 1996, 15, 1122–1125.
71. M. D. O'Donnell, S. J. Watts, R. V. Law and R. G. Hill, J. Non-Cryst. Solids, 2008, 354, 3554–3560.
72. M. D. O'Donnell, S. J. Watts, R. G. Hill and R. V. Law, J. Mater. Sci.: Mater. Med., 2009, 20, 1611–1618.
73. M. Edén, J. Non-Cryst. Solids, 2011, 357, 1595–1602.
74. W. Smith and T. R. Forester, J. Mol. Graphics, 1996, 14, 136–141.
75. I. T. Todorov, W. Smith, K. Trachenko and M. T. Dove, J. Mater. Chem., 2006, 16, 1911–1918.
76. J. D. Gale and A. L. Rohl, Mol. Simul., 2003, 29, 291–341.
77. I. Lebecq, F. Désanglois, A. Leriche and C. Follet-Houttemane, J. Biomed. Mater. Res., 2007, 83A, 156–168.
78. T. Abe, J. Am. Ceram. Soc., 1952, 35, 284–299.
79. A. C. Wright, A. G. Clare, B. Bachra, R. N. Sinclair, A. C. Hannon and B. Vessal, Trans. Am. Crystallogr. Assoc., 1991, 27, 239–254.
80. E. Ratai, M. Janssen and H. Eckert, Solid State Ionics, 1998, 105, 25–37.
81. J. D. Epping, W. Strojek and H. Eckert, Phys. Chem. Chem. Phys., 2005, 7, 2384–2389.

---

Footnote

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

---