A modified random network model for P₂O₅–Na₂O–Al₂O₃–SiO₂ glass studied by molecular dynamics simulations

Abstract

We investigated the short- and medium-range structural features of sodium aluminosilicate glasses with various P₂O₅ (0–7 mol%) content and Al/Na ratios ranging from 0.667 to 2.000 by using molecular dynamics simulations. The local environment evolution of network former cations (Si, Al, P) and the extent of clustering behavior of modifiers (Na⁺) is determined through pair distribution function (PDF), total correlation function (TDF), coordination number (CN), Q_xⁿ distribution and oxygen speciation analysis. We show that Al–O–P and Si–O–Al linkage is preferred over other connections as compared to a random model and that Si–O–Si linkage is promoted by the P₂O₅ addition, which is related to structural heterogeneity and generates well-separated silicon-rich and aluminum–phosphorus-rich regions. Meanwhile, due to the relatively high propensity of Al to both Si and P, heterogeneity can be partly overcome with high Al content. A small amount of Si–O–P linkages have been detected at the interface of separated regions. Clustering of Na⁺ is also observed and intensified with the addition of P₂O₅. Based on the simulated structural information, a modified random network model for P₂O₅-bearing sodium aluminosilicate glass has been proposed, which could be useful to optimize the mobility of sodium ions and design novel functional glass compositions.

---

1 Introduction

Chemical strengthening, achieved by ion exchange through immersing an alkali-containing glass into a molten salt bath, generates high compressive stress in thin or irregular shaped glass objects without measurable optical distortion, making it one of the most important techniques of glass strengthening methods.¹ Sodium aluminosilicate glasses (NAS) are suitable for ion exchange strengthening and currently receive significant interest not only because of their excellent mechanical properties and chemical durability, leading to wide applications in touch-screen displays,² windshields,^3,4 and architectures,^5,6 but also because they can serve as a good model system to introduce a wide range of oxides for tailoring their properties.^7–11 Among these oxide additions, P₂O₅ has been recently studied and found to be able to increase the hardness, compressive stress and depth of ion-exchange layer by causing depolymerized silicate anions.¹²

Phosphorus pentoxide (P₂O₅), even with a very small addition, is known to play an important role in aluminosilicate glass system and can impact the physical and chemical properties.^13,14 The immiscibility of phosphate-rich and silicate-rich phases caused by P₂O₅ addition is also observed,^15,16 accounting for crystallization after further thermal treatment. The ³¹P MAS NMR study of the P₂O₅ containing NAS glasses with a molar ratio of Al/Na ranging from 0 to about 1.27 shows that in peralkaline region P₂O₅ is partly present as phosphate tetrahedra attached to the aluminosilicate network charge balanced by sodium, while in peraluminous region it integrated into the glass former network as monomeric tetrahedra for which all three available single bonded oxygens are connected to tetrahedrally coordinated Al.¹⁷ In lithium aluminosilicate system, the P, Al and Si local environment investigations also show that P is mainly presented as orthophosphate and pyrophosphate species with a low amount of Al₂O₃ and gradually enters the Si and Al based glass forming network with increasing molar ratio of Al₂O₃/Li₂O.¹⁸ While in metaluminous sodium aluminosilicate glasses, P is found to mainly integrate into the glass network via P–O–Al bonding because of the alternative charge compensation scheme for [AlO₄]⁻ and [PO₄]⁺ units, and the formation of Si–O–P bonding played only a minor role.¹⁹ Although these experimental results provide valuable information for distinguishing the topology of the glassy states with phosphorus, fully understanding the structural change aluminosilicate glasses due to P₂O₅ and associated ion exchange strengthening still remains a grand challenge.

In the past few decades, the structure models for oxide glasses have been developed. The basic concept was first established in the 1930s,²⁰ and then a modified random network (MRN) model was proposed by Greaves²¹ in mid-1980s, in which the oxide glass is described as partially covalent random networks made up of corner-shared glass former-oxide tetrahedra with glass modifiers in the inter-networking regions. Recently, a new modified model has been proposed to explain the prevention of the alumina-avoidance principle in term of aluminosilicate glass system,²² in which the network structure region composed of well-separated homonuclear framework containing only silicon–oxygen connection and heteronuclear framework containing both silicon and aluminum tetrahedra, and the inter-networking region composed of non-bridging oxygens and non-framework cations. Considering the fact that the ionic diffusion mechanisms are mostly vacancy-based rather than interstitial-based,²³ and the possibility that the inter-networking region formed due to phosphorus addition, such modified random network models are instructive to understand the ion diffusion in P₂O₅–Na₂O–Al₂O₃–SiO₂ glass system.

Molecular dynamics (MD) simulations is an effective method to study the atomic structures and structure–property relationship of glass materials,²⁴ and has been successfully applied to the investigation of structural heterogeneity in borosilicate²⁵ and oxyfluoride glasses.²⁶ In the present work, to explore the detailed structural information of P₂O₅-bearing sodium aluminosilicate glass, compositions with various P₂O₅ content and various Al/Na ratios were investigated by using MD simulations. Short- and medium-range structural features, such as pair distribution function (PDF), total correlation function (TDF), coordination number (CN), oxygen speciation, linkage between different network former cations (T–O–T connections, T = Si, Al, P in this work), network former species (Q_xⁿ), and network connectivity (NC), were analyzed. Based on the structural information from simulations, a modified random network model for P₂O₅-bearing sodium aluminosilicate glass is proposed, considering the relations between different network formers and percolated channels with sodium clustering behaviors. This structural model is then used to better understand ion-exchange efficiency acceleration of phosphorus addition in sodium aluminosilicate glass from the perspective of the glass network topology.

2 Simulation details

To fully investigate the effect of P₂O₅ on the structures of sodium aluminosilicate glass network, a wide range of compositions were considered in the present research, including the phosphorus-free and phosphorous-bearing glasses with variable Al/Na ratios, n%P₂O₅–(1 − n%)[15Na₂O–xAl₂O₃–(85 − x)SiO₂] (mol%, n = 0, 5, and x ranges from 1 to 30), labeled as NA (n = 0) and PNA (n = 5), respectively, where N, A, and P refer to Na₂O, Al₂O₃, and P₂O₅, as well as a series of compositions of variable phosphorus content with fixed Al/Na ratio, x%P₂O₅–(1 − x%)[15Na₂O–nAl₂O₃–(85 − n)SiO₂] (mol%, n = 10, 15 and 20, and x ranges from 0 to 7), labeled as LAP (n = 10), MAP (n = 15) and HAP (n = 20) for peralkaline, metaluminous, and peraluminous glass compositions, respectively, where LA, MA, and HA refer to low, medium, and high Al content, and P refer to P₂O₅. Table 1 lists specific compositions and naming of the simulated glass system. The initial glass densities were calculated by using the software Glass Engineer System.^27,28 MD simulations were performed using the DL_POLY version 2.20 package.²⁹ The modified Teter potential were employed for more efficiently calculating the glass structure and reducing the computational time, which has been expanded by Du et al. and widely used in silicate,^23,30 aluminosilicate,^11,26,31,32 and phosphosilicate^24,33,34 glass systems.

Table 1 Compositions (mol%) and densities of the simulated glass system

Series	SiO2	Al2O3	Na2O	P2O5	Al/Na	Initial density (g cm^−3)	Simulated density (g cm^−3)
NA	84	1	15	—	0.067	2.350	2.376 ± 0.007
	81	4	15	—	0.267	2.367	2.388 ± 0.006
	78	7	15	—	0.467	2.383	2.411 ± 0.001
	75	10	15	—	0.667	2.398	2.418 ± 0.006
	73	12	15	—	0.800	2.408	2.418 ± 0.005
	72	13	15	—	0.867	2.412	2.438 ± 0.005
	71	14	15	—	0.930	2.417	2.432 ± 0.007
	70	15	15	—	1.000	2.421	2.441 ± 0.013
	69	16	15	—	1.067	2.435	2.447 ± 0.010
	68	17	15	—	1.133	2.449	2.455 ± 0.002
	67	18	15	—	1.200	2.463	2.466 ± 0.018
	66	19	15	—	1.267	2.473	2.461 ± 0.008
	65	20	15	—	1.333	2.482	2.489 ± 0.002
	63	22	15	—	1.467	2.502	2.491 ± 0.008
	60	25	15	—	1.667	2.532	2.514 ± 0.021
	57	28	15	—	1.867	2.562	2.522 ± 0.002
	55	30	15	—	2.000	2.581	2.537 ± 0.012
PNA	79.8	0.95	14.25	5	0.067	2.379	2.398 ± 0.010
	76.95	3.8	14.25	5	0.267	2.393	2.402 ± 0.005
	74.1	6.65	14.25	5	0.467	2.407	2.420 ± 0.013
	71.25	9.5	14.25	5	0.667	2.419	2.423 ± 0.003
	69.35	11.4	14.25	5	0.800	2.427	2.427 ± 0.007
	68.4	12.35	14.25	5	0.867	2.431	2.424 ± 0.008
	67.45	13.3	14.25	5	0.930	2.435	2.440 ± 0.007
	66.5	14.25	14.25	5	1.000	2.437	2.437 ± 0.005
	65.55	15.2	14.25	5	1.067	2.437	2.440 ± 0.012
	64.6	16.15	14.25	5	1.133	2.437	2.444 ± 0.006
	63.65	17.1	14.25	5	1.200	2.438	2.450 ± 0.005
	62.7	18.05	14.25	5	1.267	2.445	2.455 ± 0.003
	61.75	19	14.25	5	1.333	2.454	2.462 ± 0.001
	59.85	20.9	14.25	5	1.467	2.472	2.478 ± 0.003
	57	23.75	14.25	5	1.667	2.499	2.511 ± 0.001
	54.15	26.6	14.25	5	1.867	2.509	2.504 ± 0.004
	52.25	28.5	14.25	5	2.000	2.526	2.504 ± 0.003
LAP	75	10	15	0	0.667	2.398	2.418 ± 0.006
	74.25	9.9	14.85	1	0.667	2.403	2.419 ± 0.007
	73.5	9.8	14.7	2	0.667	2.407	2.421 ± 0.005
	72.75	9.7	14.55	3	0.667	2.411	2.422 ± 0.016
	72	9.6	14.4	4	0.667	2.415	2.422 ± 0.006
	71.25	9.5	14.25	5	0.667	2.419	2.423 ± 0.003
	70.5	9.4	14.1	6	0.667	2.423	2.429 ± 0.008
	69.75	9.3	13.95	7	0.667	2.427	2.437 ± 0.007
MAP	70	15	15	0	1.000	2.421	2.441 ± 0.013
	69.3	14.85	14.85	1	1.000	2.425	2.445 ± 0.005
	68.6	14.7	14.7	2	1.000	2.428	2.440 ± 0.002
	67.9	14.55	14.55	3	1.000	2.432	2.446 ± 0.007
	67.2	14.4	14.4	4	1.000	2.436	2.439 ± 0.006
	66.5	14.25	14.25	5	1.000	2.437	2.437 ± 0.005
	65.8	14.1	14.1	6	1.000	2.438	2.439 ± 0.012
	65.1	13.95	13.95	7	1.000	2.439	2.436 ± 0.011
HAP	65	20	15	0	1.333	2.482	2.489 ± 0.002
	64.35	19.8	14.85	1	1.333	2.476	2.474 ± 0.006
	63.7	19.6	14.7	2	1.333	2.471	2.476 ± 0.004
	63.05	19.4	14.55	3	1.333	2.465	2.472 ± 0.005
	62.4	19.2	14.4	4	1.333	2.459	2.448 ± 0.002
	61.75	19	14.25	5	1.333	2.454	2.463 ± 0.001
	61.1	18.8	14.1	6	1.333	2.449	2.461 ± 0.004
	60.45	18.6	13.95	7	1.333	2.444	2.454 ± 0.004

---

During the simulated melt-and-quench procedure, canonical (NVT), microcanonical (NVE) and constant temperature and pressure (NPT) ensembles have been chosen to obtain reasonable glass structure. In each simulation, around 3000 atoms were used in the simulation box with periodic boundary conditions. Larger simulation cells with 6000 atoms and smaller simulation cells with 1000 atoms have also been simulated for the compositions with R_Al/Na = 1.000 in NA and PNA series to study the effect of system size. Related data are available in ESI (see SF 1)† and show that structural data converges for systems with 3000 and 6000 atoms per simulation cell while the smaller system with 1000 atoms showed some discrepancies. Considering both the accuracy of glass model and computational efficiency, all the results and discussion reported here are based on the simulation cell with 3000 atoms. The initial atom positions were generated by randomly putting the atoms in a cubic simulation cell with cell volume consistent to the calculated glass density. To avoid atoms being placed too close to each other during the generation of the initial configurations, constraints of minimum interatomic distance based on ionic radii were applied. After relaxing at 0 K and 0 Pa, the melting temperature reached at 6000 K through 2000 and 4000 K for fully eliminating the memory effects of the initial configuration. After that, a two-stage process with a step-wise cooling strategy was used with a nominal cooling rate of 1 K ps⁻¹. First, before the temperature reached down to 1500 K through 4000, 3000, 2500, 2000 K, each temperature stair adopted NVT ensemble for 200 000 steps (200 ps), followed by NVE for another 200 000 steps (200 ps). Second, NPT and NVE runs were performed in sequence at the temperature ranging from 1500 K to 300 K through 1000, 750, 700, 600, 500, 400 K, aiming to shrink the voids in glass samples as well as relax the system to equilibrium volume. The final trajectories of glass system, recorded every 0.05 ps, were generated by additional 200 ps NVE run at 300 K. These 4000 configurations were used for final structure analysis. For each composition, three parallel simulations were conducted starting from different initial random structures and the statistics analyses were averaged among the three samples. The simulated densities are shown in Table 1 and ranges from 2.37 to 2.54 g cm⁻³, consistent with previous reported experimental result.^31,35,36 Compared with initial densities, the deviations are less than 2% and 1% in NA and PNA series, respectively, which is quiet reasonable for MD simulations.³¹

3 Results

3.1 General glass network structures

The analysis of network former species (Q_xⁿ) has been carried out as shown in Fig. 1, which is an important method to understand the intermediate range structural information in oxide glass, where n is the number of bridging oxygen (BO) coordinating with a glass network former atom x (x= Si, Al, P in the present work). [SiO₄] tetrahedra in both NA and PNA compositions mainly exist as Q_Si⁴ species (Fig. 1a), with the coexistence of Q_Si² and Q_Si³ groups, and the amount of Q_Si⁴ increasing with consuming of Q_Si² and Q_Si³ as Al/Na ratio increasing, consistent with the experimental and theoretical results from the earlier literatures,^32,37 indicating a polymerized silicon–oxygen network,^9,18,33,38 which could be attributed to more sodium ions acting as charge compensator for [AlO₄] tetrahedra. It is worthy to note that the P₂O₅ addition increases Q_Si⁴ amount with reducing Q_Si³ in peralkaline compositions (Fig. 1a). On the contrary, the Al environment has not been significantly influenced by Al content or P₂O₅ addition (Fig. 1b), with Q_Al⁴ maintaining the dominant species coexistence with a minor amount of five-fold coordinated sites Q_Al⁵.³⁹ As shown in Fig. 1c, Q_Pⁿ exhibits more complicated effects compared with Si and Al. In peralkaline region, P is mainly located as depolymerized units Q_P¹ and Q_P². Due to the fact that these units need positive charge to maintain local neutrality, sodium ions have been partly scavenged from silicon–oxygen network, leading to the increased amount of Q_Si⁴, which means a denser silicon–oxygen network in peralkaline region (see Fig. 1a). These depolymerized units disappear and Q_P³ species become dominant as increasing Al/Na ratio. After Al/Na > 1, P mainly exists as Q_P⁴ with consuming Q_P³, providing charge balance for [AlO₄] tetrahedra,¹⁹ thus a fraction of sodium ions transform the role of charge compensator to network modifier and move back to the vicinity of [SiO₄] units, resulting in a slight reduction of Q_Si⁴ and increase of Q_Si³.

Fig. 1 Q_xⁿ distribution for (a) Si, (b) Al, (c) P and (d) Q̄_xⁿ and network connectivity (NC) in the compositions with different Al/Na ratio.

To further investigate the glass network in NA and PNA, the mean Q_xⁿ was yield by

where x_n is the percentage of this type Qⁿ for the element x. Then the glass network connectivity (NC) was calculated bywhere m_x is the atom number of the element x in the simulation. As shown in Fig. 1d, all average number of bridging oxygen exhibit a similar trend with the dominant species, Q_Si⁴, Q_Al⁴, and Q_P⁴, respectively, which is an indication that the polymerization of Si–O network as well as P–O network increase with Al/Na ratio, especially for peralkaline compositions. The addition of P₂O₅ enhances the polymerization of Si–O network in terms of the higher value of Q̄_Siⁿ. However, the network connectivity of the overall glass structure exhibits a depolymerized behavior after P₂O₅ addition, even it increases with Al/Na ratio, implying a structural heterogeneity containing separated denser Si–O network with higher polymerization than overall glass structure and depolymerized P-rich regions.

The polymerized effect of P–O on Si–O network is more significant for peralkaline compositions (see SF 2† for LAP), with increasing P₂O₅ content, the amount of Q_Si⁴ increases, causing higher value of Q̄_Siⁿ. The influence gradually diminishes with increasing Al/Na ratio, that is for metaluminous and peraluminous compositions with various P₂O₅ content (see SF 3&4† respective for MAP and HAP), the decreased value of Q_Si⁴ and Q̄_Siⁿ is attributed to the partly release of sodium ions from [PO₄] units to [SiO₄] units.

In aluminosilicate glasses, phosphorus can play two different roles: first, for peralkaline compositions, P mainly exists as depolymerized species Q_P¹ and Q_P² and partly extracts sodium ions from silicon–oxygen network to maintain local charge neutrality, and second, for metaluminous and peraluminous compositions, Q_P⁴ is the dominant species with the breakage of double bond to provide an alternative charge compensation for [AlO₄] units,^14,19 which can be verified by oxygen speciation analysis as shown in Fig. 2. In peralkaline region, with increasing Al/Na ratio, the amount of non-bridging oxygen (NBO) decreases while that of BO increases (Fig. 2a), indicating a higher network connectivity. When Al/Na > 1, three-bonded oxygen (TBO) appears and increases dramatically with consuming of BO, due to the lack of sodium ions to charge compensate [AlO₄] tetrahedra. Contrary to the effect of Al, the P₂O₅ addition has increased the amount of NBO and decreased that of TBO. For BO, the effect of P₂O₅ is more complicated, that is, in peralkaline region, the P₂O₅ addition reduces the amount of BO while increases it in peraluminous region. Such trend can be explained that when Al/Na < 1 there is enough sodium ions acting as charge compensator for [AlO₄] units, with excess Na⁺ to generate NBOs in [SiO₄] units. The added P would attract Na⁺ around [SiO₄] and [AlO₄] units to provide charge balance for Q_P¹ and Q_P², thus more sodium ions act as modifiers by breaking up the glass network with increased NBO and decreased BO. When Al/Na > 1, Q_P⁴ units partly replace the role of TBO and connect with [AlO₄] units, transforming TBOs to BOs thus increasing the amount of BO with consuming TBO. These are consistent with simulated results of the LAP, MAP and HAP compositions in Fig. 2b, wherein the amount of NBO monotonically increases and TBO decreases with increasing P₂O₅ content. For LAP compositions, the amount of BO decreases to a large extent with added P₂O₅, while for MAP and HAP, the effect gradually diminishes, the amount of BO even keeps almost unchanged for HAP, attributing to the competition between the depolymerization and charge compensation effects of phosphorus. It is the same reason that the peak shape of BO in Fig. 2a becomes asymmetric and peak position shifts from the ratio of Al/Na = 1 to 1.2, at which point there is a balance between the two different effects of P₂O₅ on glass network.

Fig. 2 (a) Al/Na ratio dependence of oxygen speciation in NA and PNA, and (b) P₂O₅ content dependence of oxygen speciation in LAP, MAP and HAP.

3.2 Local environment around sodium ions

Focusing on the effect of P₂O₅ upon sodium ions, the total correlation function of Na–O for LAP compositions, as shown in Fig. 3a, assesses the local environmental evolution with added P₂O₅. The overall shape and peak position change little with composition. From the enlarged image of first peak in the inset, there is only a slight left shift of peak position, which is almost negligible. The result is nearly the same as it is in MAP and HAP compositions (see SF 5†). Based on the further analysis of TDF in Fig. 3b, Na–NBO bond length is about 2.34 Å while that of Na–BO being 2.55 Å. This suggests that the Na–O bond length becomes shorter when Na ions act as network modifiers than when they act as charge-compensators,^25,40 because of the fact that Na ions are more strongly bonded to the created NBO than to the surrounding BO when they act as charge compensators.

Fig. 3 Total correlation function of (a) Na–O, (b) Na–NBO and Na–BO for LAP glass compositions.

The coordination number (CN) of Na–O was calculated from the integral to the first peak of TDF and shown in Fig. 4. The cutoff value, determined from the first minimum distance of pair correlation function, was 3.25 Å for Na–O pairs.^41,42 The Na–O coordination number increases in the sequence of LAP, MAP, and HAP, as well as BO and TBO, wherein NBO decreases. Contrary to the effect of Al, the addition of P₂O₅ increases the amount of NBO around sodium ions and decreases that of BO and TBO, which lead into the result that total CN of oxygen around sodium ions does not change significantly. This is attributed to the different mechanism of Al and P to attract Na⁺. Due to the fact that [AlO₄] units are negative charged, sodium ions are attracted to provide charge compensation and bonded with bridging oxygen with larger bond length, causing increased oxygen coordination number. As for [PO₄] units, sodium ions are extracted from the silicon–aluminum network to form Q_P¹ and Q_P² and act as modifiers to bond with non-bridging oxygen of these depolymerized units. As the number of NBO around sodium ions is much lower than that of BO, the average local oxygen environment is dominated by BO and change little with P₂O₅ content.

Fig. 4 Coordination number of (a) Na–O, (b) Na–NBO, (c) Na–BO, and (d) Na–TBO for LAP, MAP and HAP glass compositions.

3.3 The distribution of sodium ions and correlation with network formers

In order to investigate the spatial organization of sodium ions, Na–Na PDFs in different compositions have been calculated and shown in Fig. 5. In all the pair distribution functions, we observe the existence of a well-defined peak, suggesting the existence of some spatial correlations among sodium ions, which is indicative of a trend to cluster within modifier channels and has been identified intrinsic in silicate glasses.^21,30,43 For peralkaline compositions (Fig. 5a, LAP), the first peak position of the distribution function left shifts slightly and the intensity increases significantly, indicating that the propensity for clustering increases with the addition of P₂O₅. Compared the Na–Na PDFs of different Al/Na ratios, we also conclude that the degree of clustering seems to decrease with increasing Al content, from the decreased difference of intensity or peak position in MAP and HAP glasses in Fig. 5b and c. It has been further verified by the PDF of Na–Na for P0 composition with various Al/Na ratios in Fig. 5d as the intensity of first peak monotonously decreases with increased Al content. These results suggest that the spatial distribution of Na⁺ is closely related to the interaction with different glass network former species.

Fig. 5 Pair distribution function of Na–Na for (a) LAP, (b) MAP, (c) HAP glass compositions and (d) P0 composition with various Al/Na ratios.

To further study the correlation between network former cations and sodium ions, the average number of Na⁺ around different network former cations through the integral of TDFs of Si–Na, Al–Na, and P–Na (see SF 6†), has been calculated and shown in Fig. 6. In all simulated compositions, the average number of Na⁺ around different network former cations decrease with increasing P₂O₅ content. In peralkaline region, the preference of sodium ions around network former cations is found to decrease in the sequence of phosphorus, aluminum, and silicon from Fig. 6a. Sodium ions tend to be distributed in the vicinity of [PO₄] units and act as modifiers to generate NBO,⁴⁴ followed by that of [AlO₄] and [SiO₄] units. The increasing Al content (LAP, MAP, and HAP) mainly affects the distribution of sodium ions around phosphorus and aluminum cations, as the average number Na⁺ around [PO₄] and [AlO₄] units decreased while that around [SiO₄] units remain almost unchanged in Fig. 6b. The distribution of sodium ions near Si is mainly affected by P. With increasing P₂O₅ content, sodium ions around Si decrease, indicating that sodium ions are extracted from the silicon–oxygen network, consistent with the above analysis of Q_Siⁿ distribution (see Fig. 1a).

Fig. 6 P₂O₅ content dependence of the average number of Na⁺ around different network-forming cations (a) in LAP, MAP, and HAP glass compositions and (b) for Si, Al and P with various Al content.

4 Discussion

4.1 Structural heterogeneity in phosphorus-bearing aluminosilicate glass

The simulated microstructures of LAP glasses are presented in Fig. 7. It is clearly that the glass structure is far from homogeneous. In the P₂O₅-free aluminosilicate glasses (Fig. 7a), there are two regions, a network structure region composed of well-separated aluminum–silicon and silicon-rich framework species as well as a channel region of non-bridging oxygens and network modifiers inserted into aluminum–silicon region, which is consistent with the proposed model in terms of aluminosilicate glasses.²² In the P₂O₅-bearing aluminosilicate glasses (Fig. 7b), the structural heterogeneity still exist and even more significant as the fraction of homonuclear framework predominantly containing [SiO₄] units increase after P₂O₅ addition. A fraction of P is distributed at the interface of silicon-rich and aluminum–phosphorus regions, interpreting the formation of minor P–O–Si linkage, which has been detected in previous NMR experiments.¹⁹ While for MAP and HAP compositions (see SF 7†), the effect of P on the heterogeneous distribution of network formers and modifiers seems weaker.

Fig. 7 Morphology screenshots of some typical simulation samples (a) LAP-P0, (b) LAP-P7. The screenshots show the bulk size with around 34 Å in length 6.8 Å for thickness. Color indications: red-O, blue-Si (stick model); bluegreen tetrahedra-Al, purple tetrahedra-P, yellow ball-Na (yellow shaded area is percolated channel region, P-rich area is stressed by bold dotted line in the P₂O₅-bearing glass composition).

In order to investigate the mutual affinity and avoidance between different network formers, the percentage of different T–O–T (T = Si, Al, P) connections in glass compositions with various P₂O₅ content has been calculated through detailed analysis of BO and TBO and shown in Fig. 8a. In the simulated glass system, with increasing P₂O₅, the percentage of Si–O–Si, Si–O–Al, and Al–O–Al decreases while that of Si–O–P and Al–O–P increases. The amount of P–O–P linkage is too small to be drawn in the chart. It is worthy to note that a fraction of Al–O–Al exits in our simulated structure. The result is to some extent breaking the Al avoidance principle which suggests that Si–O–Al linkages are more energetically favorable such that Al–O–Al linkages are not formed if there is a possibility for Si–O–Al to occur. However, the extent of obedience of this principle in aluminosilicate glasses is questionable.^45,46 Recently, prevention of the Al avoidance principle have been investigated by NMR spectral analysis and it has been proposed that breaching of the principle was essential for understanding the distribution of [SiO₄] and [AlO₄] units in alkaline-earth aluminosilicate glass.²² In our simulation, the quantity of Al–O–Al linkage decreases with P₂O₅ content while increase with Al content, which is consistent with previous experimental and simulated results.^32,47,48

Fig. 8 (a) Fraction and (b) degree of preferred connection of different T–O–T network linkage for LAP, MAP and HAP glass compositions various P₂O₅ content.

To quantify the propensity of different inter-tetrahedral connectivity, the degree of preferred connection, DPC_mn (where m and n mean network-former cations, Si, Al, and P for example), is defined and calculated, which is governed by

where f_mn is the ratio of the connection m–O–n to the whole network former connections. MD and RM mean the values yielded by the molecular dynamic simulation or the random model. According to the random model, the linkages between the network former cations are uniform, therefore, f_mn(RM) will be written as,where l₁ and l₂ are the network former cations, and

When DPC_mn > 1, m and n show some affinity and prior to form linkage in the network, whereas when DPC_mn ≤ 1, there are not. As shown in Fig. 8b, the distribution of DPC_mn indicates that the formation of T–O–T (T = Si, Al, P) connections is far away from random model, and shows a certain propensity. The priority of different T–O–T linkage in order from high to low is Al–O–P, Si–O–Al, Si–O–Si, Al–O–Al, and Si–O–P, revealing a strong affinity of [PO₄] tetrahedra for aluminum.^10,18 The lowest priority of Si–O–P might cause immiscibility and phase separation, that is, P could serve as nucleating agent in silicate glasses.⁴⁹ Although the percentage of Si–O–Si connection decreases with increased P content, DPC_SiSi has an increasing trend, indicating that P obviously improves the priority of Si–O–Si linkage, especially for LAP compositions, suggesting the formation of separated silicon–oxygen rich regions, which could induce the emerging of cristobalite (SiO₂) crystals after heat treating.⁵⁰ Such complex bonding geometries might arise from the different value of ionic field strength, P⁵⁺ > Si⁴⁺ > Al³⁺, leading to different NBO preference causing that [PO₄] units are mainly incorporated as ortho- and pyrophosphate units (Q⁰_P and Q_P¹) and show broad Q_Pⁿ distribution with more NBO while [SiO₄] units mainly exist as Q_Si³ and Q_Si⁴ with less NBO, and [AlO₄] units mainly exist as Q_Al⁴ and are not sensitive with composition, suggesting that Si–NBO is strongly preferred than Al–NBO, which has been verified by ¹⁷O MAS NMR experiments.⁵¹ As a result, the Si environment is significantly modified by P₂O₅ addition with an increase of Q_Si⁴ species, generating a more polymerized Si–O network. Besides, there is a strong affinity of phosphorus to bond with aluminum as tetrahedral coordinated phosphorus have a bond valence of 1.25, the bond valence of [AlO₄] species is 0.75, thus a P–O–Al connectivity leaves no charge on the bridging oxygen and local electrostatic neutrality would be satisfied.^18,52 Therefore, Si–O–Al connection has been consumed to form Si–O–Si and Al–O–P with P₂O₅ addition.¹⁶

Based on the propensity of connection between different network former cations, the addition of P₂O₅ has induced separated aluminum–phosphorus-rich and silicon-rich region, exacerbating the heterogeneity of glass structure, which might not be distinguished in emission scanning microscope (SEM).¹⁶ In contrast, due to the relatively high priority of Al–O–P and Si–O–Al, the increasing Al content, acting as the role of bridge and generating linkage like Si–O–Al–O–P, tends to generate a more homogenous glass structure,^15,39 which is supported by the decreased priority of Si–O–Si with addition of Al content (see SF 8†). Besides the heterogeneity of network former cations, clustering of sodium ions is also observed with addition of P₂O₅, especially in peralkaline region, which would have an effect on the mobility of sodium ions, as the connected Na–O polyhedra provide natural pathways for migration.^23,53

4.2 Proposal of a new structural model from MD simulations

Based on the above MD simulation results, a structural model for phosphorus-bearing sodium aluminosilicate glass has been proposed, as shown in Fig. 9. In our simulated P₂O₅-bearing glasses, the network composes of three complicated corner-sharing network formers Si, Al, and P with inhomogeneous distribution. Phosphorus, with more NBO in the first coordination shell, has a strong affinity to Al forming aluminum–phosphorus-rich regions and avoids to connect with Si, meanwhile, the clustered Na⁺ are preferred to appear in the vicinity of [PO₄] and [AlO₄] units. Therefore, Al–O–P connections can be located at the border of framework,¹⁰ leaving NBO around P to inter-connect with clustered Na⁺ to form percolated channels. Due to both relatively high propensity to P and Si, [AlO₄] units play a central role of separating silicon and phosphorus, acting as bridge, to form continuous network containing three network formers. Because of the effect of P on the promotion of Si–O–Si linkage, silicon-rich regions have been formed, and separated by percolated channel from aluminum–phosphorus-rich region. Due to the preferential interaction of sodium ions to [PO₄] and [AlO₄] units, sodium depletion layer is formed at the border of phosphorous–aluminum region and beneficial for silicon-rich region, where there is a more polymerized structure compared to the overall network connectivity. The clustered sodium ions around P and Al together with NBOs form percolated channels in the glass network. Considering the fact that sodium ions cluster around [AlO₄] and [PO₄] units and that NBOs can be mainly provided by [SiO₄] and [PO₄], there are two possible channel locations (Fig. 9): Channel location 1 is the interface of separated aluminum–phosphorus and silicon-rich region; Channel location 2 is inserted into silicon–aluminum region formed by Si–O–Al connection.

Fig. 9 A structural model proposed for P₂O₅-bearing sodium aluminosilicate glasses.

In the framework of the present topological model, it is reasonable to understand that for peralkaline aluminosilicate glasses, the phosphorus addition can significantly improve the ion-exchange efficiency since it provides new percolated channel location due to the existence and exacerbation of local heterogeneity and clustering within the network former and modifier species, thus increasing the amount of diffusion pathways to improve alkali migration. On the other hand, for metaluminous and peraluminous compositions, the effect of P diminishes gradually and does not have a significant influence on ion-exchange process, as Al improves the homogeneous distribution of alkali ions and phosphorus cations in some extent.

5 Conclusions

This MD simulation study reveals the structural heterogeneity in P₂O₅-bearing sodium aluminosilicate glass system, especially in peralkaline compositions. Due to the relatively depolymerized [PO₄] units causing Na⁺ ions being moved from silicate network, as well as the strong priority of P to Al, significant heterogeneity is favored in these compositions. Nevertheless, the propensity for such heterogeneity, including the inhomogeneous distribution of mutual connection of network formers and clustering extent of sodium ions, decreases with increasing Al content. Based on these results, a structural model has been proposed from the perspective of glass network topology to elucidate the effect of P on the acceleration of ion-exchange efficiency, which depends on the Al/Na ratio. With this model, we can further study the structure–properties relationship and design novel functional glass compositions with optimized ion mobility.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This work was supported by the 13th Five-Year National Key R&D Program (No. 2016YFB0303700), National Natural Science Foundation of China (U1809217 and 51672242), Opening Project of State Key Laboratory of Advanced Technology for Float Glass and Fundamental Research Funds for the Central Universities.

References

1. S. Karlsson, B. Jonson and C. Stalhandske, Glass Technol.: Eur. J. Glass Sci. Technol., Part A, 2010, 51(2), 41–54.
2. P. Ball, Nat. Mater., 2015, 14(5), 472.
3. A. K. Varshneya, Phys. Chem. Glasses: Eur. J. Glass Sci. Technol., Part B, 2017, 58(4), 127–132.
4. A. Vedrtnam, Proc. Inst. Mech. Eng., Part C, 2019, 233(4), 1334–1344.
5. V. I. Loganina, M. E. Sadovnikova, W. Jezierski and D. Małaszkiewicz, IOP Conf. Ser.: Mater. Sci. Eng., 2017, 245, 022023.
6. G. N. Hou, C. M. Zhang, W. B. Fu, G. S. Li, J. N. Xia and Y. X. Ping, Ceram. Int., 2019, 45(9), 11850–11855.
7. M. Engholm, L. Norin and D. Aberg, Opt. Lett., 2007, 32(22), 3352–3354.
8. S. Bhattacharyya, T. Höche, N. Hemono, M. J. Pascual and P. A. van Aken, J. Cryst. Growth, 2009, 311(18), 4350–4355.
9. J. Ding, Y. K. Chen, W. Chen, L. L. Hu and G. Boulon, Chin. Opt. Lett., 2012, 10(7), 071602.
10. S. Guo-Malloy, P. F. McMillan and W. T. Petuskey, J. Non-Cryst. Solids, 2016, 451, 77–83.
11. J. J. Zhao, X. X. Xu, P. C. Li, X. Y. Li, D. Q. Chen, X. S. Qiao, J. C. Du, G. D. Qian and X. P. Fan, J. Phys. Chem. B, 2019, 123(13), 3024–3032.
12. H. D. Zeng, L. Wang, F. Ye, B. Yang, J. D. Chen, G. R. Chen and L. Y. Sun, Front. Mater., 2016, 3, 53.
13. B. O. Mysen, Contrib. Mineral. Petrol., 1998, 133(1–2), 38–50.
14. G. D. Cody, B. Mysen, G. Saghi-Szabo and J. A. Tossell, Geochim. Cosmochim. Acta, 2001, 65(14), 2395–2411.
15. C. C. de Araujo, L. Zhang and H. Eckert, J. Mater. Chem., 2006, 16(14), 1323–1331.
16. M. Tarrago, M. Garcia-Valles, S. Martinez and D. R. Neuville, J. Environ. Manage., 2018, 220, 54–64.
17. M. J. Toplis and T. Schaller, J. Non-Cryst. Solids, 1998, 224(1), 57–68.
18. P. Glatz, M. Comte, L. Montagne, B. Doumert and L. Cormier, Phys. Chem. Chem. Phys., 2019, 21(33), 18370–18379.
19. A. Nizamutdinova, T. Uesbeck, T. Grammes, D. S. Brauer and L. van Wullen, J. Phys. Chem. B, 2020, 124(13), 2691–2701.
20. W. H. Zachariasen, J. Am. Chem. Soc., 1932, 54, 3841–3851.
21. G. N. Greaves, J. Non-Cryst. Solids, 1985, 71(1–3), 203–217.
22. A. R. Allu, A. Gaddam, S. Ganisetti, S. Balaji, R. Siegel, G. C. Mather, M. Fabian, M. J. Pascual, N. Ditaranto, W. Milius, J. Senker, D. A. Agarkov, V. V. Kharton and J. M. F. Ferreira, J. Phys. Chem. B, 2018, 122(17), 4737–4747.
23. A. N. Cormack, J. C. Du and T. R. Zeitler, Phys. Chem. Chem. Phys., 2002, 4(14), 3193–3197.
24. X. N. Lu, L. Deng, C. Huntley, M. G. Ren, P. H. Kuo, T. Thomas, J. Chen and J. C. Du, J. Phys. Chem. B, 2018, 122(9), 2564–2577.
25. M. Y. Wang, M. M. Smedskjaer, J. C. Mauro and M. Bauchy, J. Chem. Phys., 2019, 150(4), 044502.
26. J. J. Zhao, X. X. Xu, X. T. Chen, Q. Xu, Z. Luo, X. Qiao, J. C. Du, X. P. Fan and G. D. Qian, J. Eur. Ceram. Soc., 2019, 39(15), 5018–5029.
27. A. Fluegel, J. Am. Ceram. Soc., 2007, 90(8), 2622–2625.
28. F. X. Gan, Calculation of physical properties of inorganic glasses and component design (in Chinese), Shanghai Science and Technology Press, Shanghai, 1981.
29. W. Smith and T. R. Forester, J. Mol. Graphics, 1996, 14(3), 136–141.
30. J. C. Du and A. N. Cormack, J. Non-Cryst. Solids, 2004, 349, 66–79.
31. M. G. Ren, J. Y. Cheng, S. P. Jaccani, S. Kapoor, R. E. Youngman, L. P. Huang, J. C. Du and A. Goel, J. Non-Cryst. Solids, 2019, 505, 144–153.
32. Y. X. Zhao, J. C. Du, X. S. Qiao, X. Cao, C. Zhang, G. Xu, Y. Liu, S. Peng and G. R. Han, J. Non-Cryst. Solids, 2020, 527, 119734.
33. J. C. Du, L. Kokou, J. L. Rygel, Y. S. Chen, C. G. Pantano, R. Woodman and J. Belcher, J. Am. Ceram. Soc., 2011, 94(8), 2393–2401.
34. M. G. Ren, X. N. Lu, L. Deng, P. H. Kuo and J. C. Du, Phys. Chem. Chem. Phys., 2018, 20(20), 14090–14104.
35. J. Q. Wu, C. W. Lin, J. L. Liu, L. Han, H. Gui, C. Li, T. Y. Liu and A. X. Lu, J. Non-Cryst. Solids, 2019, 521, 119486.
36. T. Grammes, R. Limbach, S. Bruns, L. van Wüllen, D. de Ligny, E. I. Kamitsos, K. Durst, L. Wondraczek and D. S. Brauer, Front. Mater., 2020, 7, 115.
37. D. A. McKeown, F. L. Galeener and G. E. Brown, J. Non-Cryst. Solids, 1984, 68(2–3), 361–378.
38. R. Dupree, D. Holland and M. G. Mortuza, Nature, 1987, 328(6129), 416–417.
39. R. R. Deshpande, L. Zhang and H. Eckert, J. Mater. Chem., 2009, 19(8), 1151–1159.
40. M. Bauchy and M. Micoulaut, J. Non-Cryst. Solids, 2011, 357(14), 2530–2537.
41. Y. Xiang and J. C. Du, Chem. Mater., 2011, 23(11), 2703–2717.
42. Y. Xiang, J. C. Du, M. M. Smedskjaer and J. C. Mauro, J. Chem. Phys., 2013, 139(4), 044507.
43. C. D. Huang and A. N. Cormack, J. Chem. Phys., 1991, 95(5), 3634–3642.
44. L. Zhang and H. Eckert, J. Phys. Chem. B, 2006, 110(18), 8946–8958.
45. A. Jaworski, B. Stevensson, B. Pahari, K. Okhotnikov and M. Eden, Phys. Chem. Chem. Phys., 2012, 14(45), 15866–15878.
46. J. J. Ren, L. Zhang and H. Eckert, J. Phys. Chem. C, 2014, 118(9), 4906–4917.
47. S. K. Lee and J. F. Stebbins, Geochim. Cosmochim. Acta, 2009, 73(4), 1109–1119.
48. S. Kapoor, X. J. Guo, R. E. Youngman, C. L. Hogue, J. C. Mauro, S. J. Rzoska, M. Bockowski, L. R. Jensen and M. M. Smedskjaer, Phys. Rev. Appl., 2017, 7(5), 054011.
49. S. Kolay and P. Bhargava, J. Am. Ceram. Soc., 2019, 102(12), 7312–7328.
50. G. Wen, X. Zheng and L. Song, Acta Mater., 2007, 55(10), 3583–3591.
51. J. R. Allwardt, S. K. Lee and J. F. Stebbins, Am. Mineral., 2003, 88(7), 949–954.
52. J. J. Ren and H. Eckert, J. Phys. Chem. C, 2018, 122(48), 27620–27630.
53. B. Park and A. N. Cormack, J. Non-Cryst. Solids, 1999, 255(1), 112–121.

---

Footnote

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

---