Modelling the structural evolution of ternary phosphate glasses from melts to solid amorphous materials

The local and medium-range structural properties of phosphate-based melts and glasses have been characterized by means of first principles (density functional theory) and classical (shell-model) molecular dynamics simulations. The structure of glasses with biomedically active molecular compositions, (P2O5)0.45(CaO)x(Na2O)0.55 x (x 1⁄4 0.30, 0.35 and 0.40), have been generated using first principles molecular dynamics simulations for the full melt-and-quench procedure and the changes in the structural properties as the 3000 K melt is cooled down to room temperature have been compared extensively with those of the final glasses. The melts are characterized by a significant fraction of threefold (P3c) and fivefold (P5c) phosphorus atoms, but structural defects rapidly decrease during the cooling phase and for temperatures lower than 1800 K the system is free of underand overcoordinated species. The analysis of the structures of the glasses at 300 K shows a prevalence of the metaphosphate Q and pyrophosphate Q species, whereas the number of Q units, which constitute the three-dimensional phosphate network, significantly decreases with the increase of calcium content in the glass. The radial and angular distribution functions indicate that higher calcium concentration in the glass leads to an increase of the rigidity of the phosphate tetrahedral network, which has been explained in terms of the calcium's higher field strength compared to that of sodium. The structural characterization of the melts and glasses obtained from first principles simulations was used to assess and validate a recently developed interatomic shell-model forcefield for phosphate-based materials. For all three compositions, our potential model is in good agreement with the first principles data. In the glass network, the forcefield provides a very good description of the split between the shorter distances of phosphorus to non-bonded oxygen and the longer distances of the phosphorus to bonded oxygen; the phosphorus–phosphorus medium-range distribution; and the coordination environment around the Na and Ca glass modifiers. Moreover, the distribution of the Q species in the melts and glasses is in excellent agreement with the values extracted from the first principles simulations. In contrast, simulations using standard rigid ion potentials do not provide a satisfactory description of the local short-range structure of phosphate-based glasses and are therefore less suitable to model this class of multicomponent amorphous system.


Introduction
Biomedical materials for use in the body used to be dominated by metal implants, which are designed to have physical properties close to the replaced hard tissue and minimize the toxic response in the host, i.e. bio-inert materials without signicant interactions with the surrounding tissue. 1 However, one of the most important advances in this eld was the development of biomedical implants that play an active, rather than a passive, role in tissue regeneration and degrade aer the tissue has healed.This class of implant has been described as a "third generation" bio-active material. 2One commercially available bioactive glass is BioglassÒ, composed of 45.0 wt% SiO 2 , 24.5 wt% CaO, 24.5 wt% Na 2 O and 6.0 wt% P 2 O 5 , 3,4 which has been remarkably successful in many clinical applications, especially in dental and orthopaedic elds.However, questions relating to the long term effect of silica, 5 and the slow degradation of these glasses, taking one to two years to disappear from the body, 6 have prompted a search for more benign materials.
An alternative class of "third-generation" biomaterial are phosphate-based bioactive glasses, 7,8 containing PO 4 units as network formers and sodium and calcium ions as network modiers.These glasses have a chemical composition that is similar to the surrounding tissue and a dissolution rate that can be controlled by changing the phosphorus-calcium-sodium ratio. 9For example, experimental evidence shows an inverse relation between calcium content and dissolution rate, 10 which suggests that the interaction of the Ca ions with the glass network controls the glass degradation.As the solubility of a glass is linked to its structure, insight into the effect of calcium and sodium ions on the P 2 O 5 network is of fundamental interest to the design of phosphate bioglasses for applications in tissue engineering.
In this respect, molecular dynamics (MD) simulations allow for an atomic-level resolution of the structure and dynamics of the simulated system.This technique can thus be extremely useful in complementing and supporting experiments towards a fundamental understanding of the properties of amorphous materials.8][19][20][21][22][23][24][25] However only very recently, a polarizable interatomic forceeld has been developed to model ternary CaO-Na 2 O-P 2 O 5 systems. 26An alternative computational approach is the ab initio molecular dynamics (AIMD) technique, 27 which is usually based on density functional theory. 28First principles simulations are signicantly more computationally intensive than classical MD, and hence restricted to much smaller models (approximately 150 atoms) and simulation times (up to several tens of picoseconds). 29,30This methodology explicitly takes into account the electronic structure of the system and enables the accurate modelling of many-body and polarization effects.Therefore, AIMD represents a better computational approach to investigate challenging systems, such as ternary and quaternary bioglasses. 11AIMD methods have indeed been used to model the structural, electronic and dynamical properties of bioactive silicate [31][32][33][34] and, more recently, binary, ternary and quaternary phosphate glasses. 35 common approach in AIMD simulations of glasses is to conduct the quench-from-the-melt using classical MD, and the glass structure thus generated is then used as the starting point for an AIMD run.20,33,34,36 This protocol allows a partial relaxation of the structure at quantum mechanical level, but other structural features, such as the inter-tetrahedral connectivity and ionic coordination shells, depend on the quality of the force-elds used to generate the glass.Alternatively, the glass structure can be obtained using AIMD for the whole melt-andquench procedure.This full ab initio procedure is computationally very demanding and much less common than the mixed classical/ab initio approaches.32,35,37 However, the application of the AIMD melt-and-quench procedure allows an accurate investigation of the mechanism of glass formation, and the full ab initio glass produced using this approach is characterised by a quantum mechanical accuracy that extends beyond the short range.Therefore, the structural evolution of ternary phosphate systems, as a function of temperature obtained from ab initio simulations, can be used to assess the quality of interatomic forceelds for modelling the melts and solid amorphous states of this important class of materials.However, it is also important to note that the small simulation boxes that are accessible by AIMD prevent a good statistical analysis of the local-and medium-range structural properties.Using a combination of AIMD simulations and NMR experiments for uorinated bio-glasses, Pedone and co-workers have shown that the cooling rates used in full ab initio melt-andquench procedures can lead to glasses with computational defects.38 This disadvantage of AIMD simulations was one of the main reasons why we have developed an accurate shell model potential for PBGs and carried out the comparative study between the two methods reported here.
In this work, we present AIMD simulations of three ternary phosphate-based systems with biomedically applicative molar compositions, (P 2 O 5 ) 0.45 (CaO) x (Na 2 O) 0.55Àx (x ¼ 0.30, 0.35 and 0.40), both in the melt and in the amorphous state.The study focuses on the structural and dynamical differences between the liquid and the corresponding glass, on the changes of these properties during the glass formation mechanism, and on the effect of the calciumsodium ratio to the glass structure.The ab initio melt and glass structures have then been used to assess and validate the formal charge shell model interatomic potential that we have recently developed to simulate phosphate-based materials. 26

Computational methodology
Ab initio molecular dynamics AIMD simulations were carried out using the Car-Parrinello code included in the Quantum-ESPRESSO package, version 4.0.1. 39The electronic structure was described by the Perdew-Burke-Ernzerhof (PBE) gradient-corrected density functional. 40Vanderbilt ultraso pseudopotentials 41,42 represented core-valence electron interactions for all atomic species; semicore shells were explicitly included for Na and Ca.The electronic wavefunctions were expanded in a plane wave (PW) basis set with a kinetic energy cutoff of 40 Ry and k-sampling was restricted to the G point of the Brillouin zone.The time step for simulations was set to 0.121 fs and the electronic mass was set to 700 a.u.The temperature was controlled using a Nosé-Hoover chain thermostat. 43Previous ab initio simulations 20,33,44,45 have proven the accuracy of the method employed in this work in determining structural, dynamical, and electronic properties of glasses.

Classical molecular dynamics
Classical molecular dynamics (MD) simulations were carried out using the DL_POLY code (version 2.20). 46The Evans thermostat 47 was used and the timestep between successive integration of the Newtonian equations of motion was set to 2 fs.We have used a full charge shell model (SM) interatomic potential that was recently parameterized 26 to reproduce experimental structural 48 and rst principles (density functional theory) mechanical properties of crystalline phosphorus pentoxide [o 0 (P 2 O 5 ) N ]. 49Short-range Buckingham potentials were used to describe the interaction of the oxygen shells with each other and with P, Na and Ca.Three-body harmonic potentials were used to control the inter-tetrahedral O-P-O and P-O-P angles.The potential terms and parameters used in this work are reported in Table 1.We have also used the standard rigid ion (RI) partial charge potential of Teter for comparison, with short-range Buckingham terms for the pair interactions between oxygen and P, O, 19 Na and Ca. 50

Simulation protocol
The phosphate glasses considered in the present study have the chemical formula (P 2 O 5 ) 0.45 (CaO) x (Na 2 O) 0.55Àx , with x ¼ 0.30, 0.35 and 0.40 (P45C30N25, P45C35N20 and P45C40N15, where P 2 O 5 ¼ P, CaO ¼ C and Na 2 O ¼ N).The size and composition of the simulated glass systems are reported in Table 2. Previous studies have shown that the local structures and dynamics of silicate glasses can be reproduced accurately using a periodic cell of $10 Å. 36,51 The compositions that we have investigated in this work are therefore built for a similar cell, according to the experimental densities at room temperature. 52he structures of the phosphate glasses P45C30N25, P45C35N20 and P45C40N15 were generated using a full AIMD melt-and-quench procedure.For each glass, the initial conguration was created by randomly placing the atoms in the simulation box, taking care to avoid signicant overlaps by means of appropriate distance cut-offs.This conguration was relaxed in the NVE (constant number of particles N, constant volume V, and constant energy E) ensemble and, aer an initial quench of all ion velocities, the system was spontaneously (i.e.without external thermostatting) heated up to 1500 K.At this point, a series of 10 subsequent NVT (constant number of particles, volume and temperature) runs of 10-15 ps each were conducted, whose target temperatures were set at 300 K intervals from 3000 K down to 600 K, with the coefficients of the electronic wavefunction re-optimised at each new temperature.This cooling phase corresponds to a nominal cooling rate of 20 K ps À1 , which is higher than typical cooling rates accessible in classical MD simulation, but comparable or lower than those used previously in the generation of silicate 51 or phosphosilicate 32,37 glasses, using a full ab initio melt and quench approach.Finally, at 300 K the production phase was carried out for 20 ps duration.The generation of glasses using classical MD employed the same protocol.

Short-range structure -AIMD
The radial distribution functions (RDFs) for the 3000 K melts and the glasses at 300 K of the (P 2 O 5 ) 0.45 (CaO) x (Na 2 O) 0.55Àx (x ¼ 30, 35 and 40) systems, simulated using AIMD, are shown in Fig. 1.The positions and intensities of the peaks are also reported in Table 3.The positions of the peaks do not change signicantly with temperature but, due to increased thermal motion and structural disorder in the melt, the intensities of the peaks in the glasses increases by a factor of 1.5-3, compared to the melts with a corresponding increase in the full-width-halfmaximum (FWHM) values.The P-O RDFs of the glasses are Three-body (the pivot is rst atom):  characterised by two well-dened peaks at 1.52 Å and 1.61-1.64Å, the separation of which increases slightly with increasing Ca/ Na ratio (see Table 3).4), in good agreement with the theoretical tetrahedral value of 109.5 .The distribution of the P-O and O-O distances (Fig. 1) and of the O-P-O angular distributions (Fig. 3) thus denotes a structural order within each tetrahedron that is essentially constant with the change in composition.In the melts at 3000 K, the ADFs display a signicantly broader prole with tails extending below 100 for compositions P45C35N20 and P45C40N15, which will be discussed in greater detail.The P-P distances and P-O-P angles of the glasses, which reect the interconnectivity between the phosphate tetrahedra, are affected by the glass composition: FWHM of the P-P distances is equal to 0.3 Å in P45C30N25 and reduces to 0.12 Å in the glass with 40% CaO content.The FWHM of the P-O-P angles decreases from the value of 24 in the P45C30N25 glass, to 16 in P45C40N15.This observation indicates that the effect of   calcium ions on the structure of the glass is to increase the rigidity of the phosphate P-O-P interconnectivity, which could simply be interpreted in terms of the more rigid coordination shell and higher eld strength of calcium, 53 compared to that of sodium. 54,55The increased rigidity of the phosphate network with the addition of calcium content could then be linked to the observed slower solubility of glasses with higher CaO content. 10owever, it is important to note this conclusion only relies on the structural analysis of the bulk and that in order to verify this hypothesis, the processes that occur at the interface between phosphate bio-glasses and water should also be considered, 31 which will be the subject of future work.In Fig. 4, the P45C40N15 O-P-O and P-O-P angles of the 3000 K melt have been decomposed into the contributions from the twofold (P 2c ), threefold (P 3c ), fourfold (P 4c ) and vefold (P 5c ) coordinated P and O atoms.The P-O-P angle is exclusively represented by the P 4c -O 2c -P 4c angle as the contribution from the angles involving threefold and vefold coordinated P atoms are not signicant.This suggests that the broadening of the P-O-P distribution in the melt, compared to that of the glass, is only related to a distortion due to an increased thermal motion, and not to specic under-or over-coordination of the atoms in the phosphate chains.A similar situation occurs for the O-P-O angular distribution, where contributions to O-P-O from O 2c -P 4c -O 2c and O 1c -P 4c -O 1c sequences are centred at 100 and 118 respectively, comparable to the structural geometry found in crystalline o 0 (P 2 O 5 ) N . 48ig. 5 reports the ADFs of the network modiers (O-Ca-O and O-Na-O).In the melts, the O-Na-O and O-Ca-O distributions display no compositional dependence with rst peak positions at approximately 57 .A typical arrangement of oxygen atoms around the calcium and sodium atoms in the glass is reported in Fig. 6.Modiers coordinate mostly with non-bridging oxygen (NBO) atoms belonging to different PO 4 tetrahedra, hence controlling the folding and connectivity of the phosphate network.The O-Na-O distributions in Fig. 5 are characterized by quite a wide range of allowed angles, which can be explained in terms of the higher exibility of the Na cation compared with Ca.For all compositions, the O-Na-O ADFs have a peak at 50-60 , which is a result of the coordination by the Na to two oxygen atoms (possibly one NBO and one BO) belonging to the same tetrahedron (O1-Na-O2 in Fig. 6).A second peak at 90 results from the coordination with NBO of different tetrahedra (O2-Na-O3 in Fig. 6).O-Ca-O second peak positions decrease from 80 in P45C30N25 to 75 in P45C40N15.Finally, the FWHM values for the second q max (degree) 108.9 123.8 p(q max ) 4.4 Â 10 À4 0.6 Â 10 À4 300 q (degree) 109.3 128.0 p(q max ) 9.3 Â 10 À4 1.2 Â 10 À4 P45C35N20 3000 q (degree) 107.5 123.4 p(q max ) 4.4 Â 10 À4 0.6 Â 10 À4 300 q (degree) 108.0 132.4 p(q max ) 11.0 Â 10 À4 1.7 Â 10 À4 P45C40N15 3000 q (degree) 108.4 129.3 p(q max ) 6.9 Â 10 À4 1.1 Â 10 À4 300 q (degree) 109.0 127.1 p(q max ) 9.5 Â 10 À4 1.8 Â 10 À4 O-Ca-O ADF peaks are 29 and 21 for the P45C30N25 and P45C40N15 compositions respectively, suggesting that the increase of calcium content in the PBG decreases the exibility of the coordination shell of Ca and consequently increases the rigidity of the glass.
Medium-range structure -AIMD Q n distribution and network connectivity.The PO 4 tetrahedron represents the basic building block of the PBG structures.The phosphate tetrahedra are classied in terms of the number of bridging oxygens (BOs) they share with other phosphate tetrahedra in the glass. 56The various types of phosphate tetrahedra that result from this classication are labelled using the Q n terminology, where n represents the number of BOs per PO 4 tetrahedron and ranges from 0 to 3. The structure of vitreous P 2 O 5 consists of only Q 3 species that form a three-dimensional network but the inclusion of modifying metal oxides, like Ca and Na, leads to the breaking of P-BO-P bonds and their replacement by P-NBO/M (where NBO stands for nonbridging oxygen and M is the modier ion, here calcium or sodium).
The distributions of the Q n species for the three compositions are reported in Table 5 and displayed in Fig. 7.In the P 2 O 5 -CaO-Na 2 O melts and glasses, the coordination statistics denote a shi of the Q n distribution from 100% Q 3 in vitreous P 2 O 5 to lower values of n in the simulated (P 2 O 5 ) 0.45 (CaO) x (Na 2 O) 0.55Àx (x ¼ 0.30, 0.35 and 0.40) systems.Compared with the Q n distribution of the glasses, the melts display a decrease in the fraction of Q 2 species and the fraction of Q 1 species is higher (see Table 5).For all three glasses, the calculated Q n distributions show a prevalence of metaphosphate Q 2 , forming chain-like phosphate fragments and pyrophosphate Q 1 species, in agreement with the experimental 31 P NMR analysis of PBGs with similar composition. 9Table 5 also indicates that the percentage of Q 3 species is considerably smaller than Q 1 and Q 2 : Q 3 units are approximately 11% of the total PO 4 tetrahedra in the P45C30N25 glass and this value reduces to zero in P45C35N20 and to only 5.6% in P45C40N15.It is interesting to note that the glass with the intermediate composition, P45C35N20, is characterized by the highest concentration of Q 2 species, which appears to occur at the expense of the Q 3 species.Classical MD simulations of the binary NaPO 3 glass found a distribution of Q 3 units of 25.4%, which is considerably larger than the percentage found in the ternary glasses.This would suggest that the calcium ion generally has a stronger tendency to break P-BO-P bonds and form P-NBO/M species than sodium, and these modify the three-dimensional Q 3 network into a chain-like phosphate network, as illustrated by the snapshot of the AIMD simulation of the P45C30N25 glass in Fig. 6.The population of isolated Q 0 species, also referred to as orthophosphate units, decreases from 5.6% in P45C30N25 to zero orthophosphate species in the PBGs with 35% and 40% content of CaO (see Fig. 7).It is known from experiment that glasses with higher calcium content show a slower solubility than glasses with higher sodium content, 10 whereas the results from our simulations suggest that the concentration of isolated PO 4 units in the glass decreases with calcium content, which would suggest that the solubility trends of PBG could not be structurally dominated by the presence of orthophosphate species.However, it is worthy of note that due to the small simulated system sizes, absolute Q n abundances should not be considered representative and medium-range structural trends are subject to large uncertainties.
Fig. 8(a) reports the changes in the Q n distribution as the 3000 K melt of the P45C40N15 system is gradually cooled down to 300 K, and shows a decrease in Q 0 , Q 1 and Q 3 with a corresponding increase in Q 2 .Results from classical MD simulations are also presented and comparative analysis will follow in later sections.
The coordination statistics reported in Table 5 show that, in the melts, the concentration of under-and over-coordinated species is signicant (10.1-15.8%for P 3c and 1.9-3.5% for P 5c ), whereas in the glasses an ideal tetrahedral (P 4c ) environment for phosphorus is reported.Finally, Fig. 9 reports the change in P and O coordination as the melt is gradually cooled from 3000 K to 300 K, which shows that no defects are present in the system below 1500 K and that the P atoms have the tendency to be under-coordinated in the melt, rather than over-coordinated.
Coordination environment.The cumulative coordination numbers of P-P, O-O, Na-O, Ca-O and P-O neighbours for the 3000 K melt and glass at 300 K of the P45C40N15 composition are displayed in Fig. 10.The variation of the P-O, Na-O and CaO coordination number with temperature is reported in Table 6.The liquid-like character of the melt is evident from the relatively continuous trend of the n(r) curves.In particular, the stepwise character of the P-P running coordination number in the glass indicates a perfect P tetrahedral environment.Each P atom is linked (through bridging oxygens) to two other P atoms in the glass, whereas the melt displays a more disordered intertetrahedral arrangement, where the P atoms have higher probability to be in a P coordination environment that differs from two.The Na-O and Ca-O pairs also indicate a higher degree of disorder in the melts, as their running coordination numbers monotonically increase from the glass to the melts (Fig. 10).Table 5 Coordination statistics of the phosphorus network from the AIMD 3000 K melts and glasses at 300 K, calculated using a 2.3 Å cutoff to define the coordination shell of P: Q n distribution, network connectivity (NC) and phosphorus (P) coordination.AIMD simulations

Composition
T (K) However, the liquid-like character of the melt does not emerge from the running coordination numbers of the P-O pair, which are still well dened in the melt and not signicantly different from the glass.However, a detailed analysis of the total fraction of under-and over-coordinated P atoms (see Fig. 9) shows that, between 3000 K and 1800 K, there are a signicant percentage (3-11%) of 3-coordinated P species.The prole of the O-O coordination number is also not signicantly affected by the melting process (see Fig. 10).The analysis of the under-and over-coordinated species shows that the oxygen atoms are mostly 2-coordinated, with less than 2% of O 0c species and no 3coordinated oxygen atoms at 3000 K (see Fig. 9).Table 7 summarizes the average oxygen coordination shells of P, Na and Ca for the glasses at 300 K and their decomposition into the BO and NBO components.

Paper
Journal of Materials Chemistry B ratio for the decomposed P-O shows a slightly higher NBO fraction, which reects the prevalence of meta-phosphate.The decomposition of the Na-O and Ca-O into BO and NBO components in Fig. 11 shows that the sodium and calcium are mostly coordinated to the non-bonding oxygens (see Table 7) and that the M-BO and M-NBO distances have two quite distinct values: the Na-NBO and Ca-NBO pairs have well dened peaks at 2.32 Å and 2.35 Å, respectively; Na-BO has a broad peak approximately centered at 2.6 Å; the Ca-BO distance is between 2.6 Å and 2.8 Å in the glasses with composition P45C35N20 and P45C40N15.The P45C30N25 Ca-BO RDF does not have any rst coordination shell peak, indicating that for this glass composition, the calcium is not coordinated to bonding oxygen.

Comparison of AIMD and classical MD melts and glasses
In this section we report a comparison of the structure of the phosphate melts and glasses obtained from the ab initio and    classical MD melt-and-quench procedures.The structural information obtained using rst principles simulations can be used as a direct reference to assess and validate the shell model (SM) interatomic forceeld.
The RDFs and O-P-O ADF of the glass with composition P45C35N20 are shown in Fig. 12.The SM provides general good agreement with the ab initio data.In particular, the forceeld is able to simulate the split between the P]O and P-O and gives a good description of the intra-tetrahedral geometry, as evidenced by the O-P-O angle distribution in Fig. 12.The most signicant difference is found in the P-P distribution, which in the SM glass has a signicantly more intense rst peak than the ab initio glass.This is most likely due to the strong force constant in the threebody P-O-P potential (see Table 1), which was derived to obtain a good description of the inter-tetrahedral P-O-P angle and mechanical properties of crystalline o 0 (P 2 O 5 ) N . 49 comparison of the P-O, Na-O and Ca-O distances in the glasses (P 2 O 5 ) 0.45 (CaO) x (Na 2 O) 0.55Àx (x ¼ 0.30, 0.35 and 0.40), decomposed into the BO and NBO components, is reported in Table 8, together with the short-range structural data for P45C30N25 obtained from XRD experiments. 52Also reported are the structures of the glasses simulated with the rigid ion (RI) potential model, using a modied Teter potential.However, for the systems with composition P45C30N25 and P45C40N15, it was not possible to generate the structures of the glasses at 300 K from a full melt-and-quench approach; during the high temperature molecular dynamics simulations, the temperature of these systems did not equilibrate to the target temperature but their energies and temperatures diverged to very large unphysical values.Therefore, for P45C30N25 and P45C40N15, we have used the glass structure generated with the shell model potential to carry out 40 ps of NVT (300 K) simulations, using the RI model from which the P-O, Na-O and Ca-O distances were extracted.Table 8 shows that the split between P-BO and P-NBO in the SM glasses (0.09-0.15 Å) is similar to that in the AI glasses (0.12-0.13 Å).This classical potential predicts an average P-O distance, which is slightly shorter, but very close to the experimental estimates (see Table 8).We can therefore conclude that the SM potential provides an  8) but predicts a difference between P-BO and P-NBO distances of 0.06 Å, approximately one-half of that found in the SM and AI glasses and experiment.The mixed character of the P-O bonds, whose ionic-covalent balance depends on the local coordination environment, such as the composition and geometry of the Na/Ca cations surrounding the PO 4 tetrahedra, 11 implies that a satisfactory description of the localand medium-range structural features of ternary phosphate-based glasses requires the application of a shell model approach to explicitly include the oxide polarization in the potential.
The coordination environment of the modiers Na and Ca is also described fairly accurately by the SM (see Table 8).In particular, for the sodium ion, the SM predicts a coordination shell of 7.0-7.4,slightly larger than in the ab initio glasses (6.4-6.5),whereas for the calcium ion the agreement is much closer: 6.8-6.9 from SM calculations, compared with 7.2 from AI simulations.Moreover, the SM correctly and quantitatively predicts that Na and Ca ions coordinate mostly to non-bonding Table 9 Coordination statistics of the phosphorus network for the 3000 K melts and glasses at 300 K obtained from AI and SM molecular dynamics simulations: Q n distribution, network connectivity (NC) and phosphorus (P) coordination As previously stated, Fig. 8 shows the Q n distribution evolution as a function of time for each methodology, using the same initial congurations and analogous simulation protocols.Results for both AIMD and classical SMMD methodologies show good agreement.For P45C40N15, the Q 2 proportion from AIMD results increases from 40.1% in the 3000 K melt to 66.7% at 300 K, compared with 46.3% to 66.7% respectively using classical MD.Similar good agreement is shown for the temperature trends and absolute values of all other Q n species.
For simulations of bioactive phosphosilicate glasses and melts, Tilocca 58 has shown that SM potentials also perform better than RI models in these systems.The correct description of the dynamical balance between the interconversion of Q n species during the cooling of a melt is achieved with SM as evaluated by higher level reference data. 58inally, Table 9 reports a comparison of the distribution of the Q n species for the three compositions, both in the melt and in the glass, computed using the AI and SM models.The small samples considered in the present study clearly do not provide an exhaustive test of the accuracy of the SM potential to describe Q n abundances; however, they still provide an adequate representation of the local dynamics of the individual Q n species. 58Table 9 shows that for both melts and glasses, and for all three different compositions, the distribution obtained using the polarizable forceeld follows very closely the one obtained from ab initio calculations.These results therefore indicate that the full charge shell model interatomic potential that we have developed and used in this study is capable of modelling the structure of ternary phosphate systems both in the melt and in the amorphous state.

Conclusions
The local and medium-range structures of ternary phosphatebased systems with composition (P 2 O 5 ) 0.45 (CaO) x (Na 2 O) 0.55Àx (x ¼ 0.30, 0.35 and 0.40), both in the melt and in the amorphous state, have been investigated using ab initio and classical molecular dynamics simulations.The structure of the glasses has been obtained using a full melt-and-quench procedure.The changes in the structural properties, as the 3000 K melt is cooled down at room temperature, have been compared extensively with that of the nal glasses.The melts are characterized by a signicant fraction of threefold (P 3c ) and vefold (P 5c ) phosphorus atoms but structural defects rapidly decrease during the cooling phase.At temperatures lower than 1800 K, the system is free of under-and over-coordinated species.Structural analysis of the glasses at 300 K shows a prevalence of metaphosphate Q 2 and pyrophosphate Q 1 species, whereas the number of Q 3 units, which constitute the three-dimensional phosphate network, signicantly decreases with increased calcium content.Calculation of the pair and angular distribution functions suggests that the rigidity of the phosphate tetrahedral glass network increases with the concentration of calcium, an observation which is interpreted in terms of the tendency of the calcium ion to be a stronger coordinator than sodium.
Moreover, because the AIMD quench from the melt allowed the full relaxation of the medium-range structure at the ab initio level, it provides us with a very accurate reference structure of the glass, which has been used to assess the quality of a formal charge shell model interatomic potential.The former has been parameterized from experimental and density functional theory calculations of structural and mechanical properties of crystalline o 0 (P 2 O 5 ) N .For all three compositions, our potential model gives good agreement with the ab initio data.In particular, the shell model simulates a value for the splitting between the P-BO and P-NBO distances (0.09-0.15 Å) very close to that obtained for the AI glasses (0.12-0.13 Å); the SM and AI glasses have similar intra-tetrahedral O-P-O angle distributions; the coordination shell of M-BO and M-NBO (M ¼ Na or Ca) are in close agreement.On the other hand, the use of a standard rigid ion potential is unable to provide a satisfactorily accurate description of the local short-range structure and, as previously shown by Tilocca for soda-lime phosphosilicate systems, 58 the inclusion of polarization in the interatomic potential is necessary in order to model multicomponent glasses.

Fig. 2
reports the temperature evolution of the P-O RDFs for P45C40N15, where it is noted that the P-O splitting occurs at 600 K.The angular distribution functions (ADFs) of the tetrahedral network (O-P-O and P-O-P) are displayed in Fig.3and the positions and intensities of the intra-tetrahedral (O-P-O) and inter-tetrahedral (P-O-P) ADF peaks are reported in Table 4.In the glasses, the O-P-O angle distributions are quite narrow and, depending on the composition, centred at 108.0-109.3degrees (see Table

Fig. 2
Fig. 2 Radial distribution functions of the P-O atomic pair for P45C40N15 at different temperatures.AIMD simulations.

Fig. 4
Fig. 4 Decomposition of the O-P-O (top) and P-O-P (bottom) angular distribution functions into the contributions from the twofold, threefold, fourfold and fivefold coordinated P and O atoms of the 3000 K melt with composition P45C40N15.AIMD simulations.

Fig. 6
Fig. 6 Snapshot from the AIMD trajectory of the P45C30N25 glass showing the typical calcium and sodium coordination shells.Ca, Na, O and P are pictured as green, purple, red and pink spheres, respectively.

Fig. 7
Fig. 7 Distribution of Q n species in the 3000 K melts (top) and glasses at 300 K (bottom) computed from the AIMD simulations of the (P 2 O 5 ) 0.45 CaO x Na 2 O 0.55Àx (x ¼ 30, 35 and 40) systems.

Fig. 9
Fig. 9 Average total fraction of threefold-and fivefold-coordinated P, threefold coordinated O and uncoordinated O computed from the AIMD simulations melt of the P45C40N15 system.

Fig. 10
Fig. 10 Cumulative (integrated) coordination numbers calculated for the 3000 K melts and glasses at 300 K computed from the AIMD simulations melt of the P45C40N15 system.

Fig. 11
Fig. 11 Glass at 300 K with composition P45C30N25.Decomposition of the P-O radial distribution functions into their BO (solid line) and NBO (dashed line) components (top).The inset reports the decomposition of the P-NBO radial distribution function into their P Q 1 -NBO and P Q 2 -NBO components.Decomposition of the Na-O and Ca-O into their BO and NBO components (bottom).AIMD simulations.

Fig. 12
Fig. 12 Comparison of the Ca-O, Na-O, O-O, P-O and P-P radial distribution functions and O-P-O angular distribution functions of the P45C35N20 glass at 300 K obtained from full ab initio (AI) and shell model (SM) molecular dynamics melt-and-quench procedures.

Table 1
Potential terms and parameters used in this work min ¼ 0.0/r max ¼ 5.0

Table 3
Positions and amplitudes for 3000 K melts and glasses at 300 K from AIMD simulations Composition T (K) P-O N a -O Ca-O O-O P-P P45C30N25 3000 r max ( Å) 1.

Table 6
Coordination numbers of phosphorus, sodium, and calcium to oxygen at various temperatures of the glass with composition P45C40N15, calculated using the cutoff specified in the table.AIMD simulations

Table 7
Coordination environments of the phosphorus, sodium, and calcium to oxygen for glasses at 300 K decomposed into their BO and NBO components.Average coordination numbers (CN av ) are obtained using P-O, Na-O and Ca-O cutoffs of 2.3 Å, 3.3 Å and 3.3 Å, respectively.AIMD simulations

Table 8
Short-range structure of the glasses at 300 K obtained from ab initio (AI) molecular dynamics, and classical molecular dynamics using the shell model (SM) and rigid ion (RI) interatomic potentials.The distance of the atom pairs corresponds to the maximum of the radial distribution function.Values within round brackets are coordination numbers obtained using P-O, Na-O and Ca-O cutoffs of 2.3 Å, 3.3 Å and 3.3 Å, respectively.Values within square brackets are the full widths at half maximum of the radial distribution functions