PCCP View Article Online

Published on 14 July 2014. Downloaded by University of Connecticut on 28/10/2014 18:30:36.

PAPER

Cite this: Phys. Chem. Chem. Phys., 2014, 16, 21135

View Journal | View Issue

On the structure of biomedical silver-doped phosphate-based glasses from molecular dynamics simulations† Richard I. Ainsworth, Jamieson K. Christie and Nora H. de Leeuw* First-principles and classical molecular dynamics simulations of undoped and silver-doped phosphatebased glasses with 50 mol% P2O5, 0–20 mol% Ag2O, and varying amounts of Na2O and CaO have been carried out. Ag occupies a distorted local coordination with a mean Ag–O bond length of 2.5 Å and an ill-defined first coordination shell. This environment is shown to be distorted octahedral/trigonal bipyramidal. Ag–O coordination numbers of 5.42 and 5.54–5.71 are calculated for first-principles and classical methodologies respectively. A disproportionation in the medium-range phosphorus Qn distribution is

Received 7th February 2014, Accepted 14th July 2014

explicitly displayed upon silver-doping via CaO substitution, approximating 2Q2 - Q1 + Q3, but not on

DOI: 10.1039/c4cp00574k

oxygen partial pair-correlation function is strong evidence for a bulk structural mechanism associated

www.rsc.org/pccp

tion is known to decrease dissolution and we show this to be a result of Ag’s local bonding.

silver-doping via Na2O substitution. An accompanying increase in FWHM of the phosphorus to bridging with decreased dissolution rates with increased silver content. Experimentally, Ag2O 2 Na2O substitu-

I. Introduction The doping of biomedical phosphate-based glasses (PBG) with species that have known biocidal effects, such as silver and copper ions, offers an alternative method for the treatment of infection.1,2 The ions can be incorporated into the glass matrix, thus providing a controlled, site-specific delivery system predicated on the degradation characteristics of the bulk glass. The bactericidal activity of doped PBG in vitro and in vivo relies on a continuous supply of the biocidal dopant into solution.3,4 Ag-doped PBG has a broad spectrum of bactericidal activity against a range of bacteria5 and associated biofilms,3 which is linked to the effect of the Ag+ ions on the bacterial respiratory chain and a loss of bacterial membrane integrity.6 Increasing the concentration of Ag in the range (P2O5)0.50(CaO)0.30(Na2O)0.20x(Ag)x (x = 0.00–0.05), leads to a decrease in glass dissolution from 1.6 to 1.1 mg mm2 h1 across the compositional limits.5 Similarly, for compositions x = 0.10, 0.15 and 0.20, dissolution rates of 1.22, 0.41 and 0.42 mg mm2 h1 were obtained respectively.3 Silver ion release rate is correlated to these degradation rates (0.076, 0.054 and 0.056 ppm h1 respectively) corresponding to a reduction in bactericidal activity with increasing Ag+ content. The dissolution behaviour and associated ion release are likely

Department of Chemistry, University College London, 20 Gordon Street, London, WC1H 0AJ, UK. E-mail: [email protected] † Electronic supplementary information (ESI) available. See DOI: 10.1039/ c4cp00574k

This journal is © the Owner Societies 2014

governed by a complex set of coupled properties such as bulk and surface structure, ionic diffusion and dissolution kinetics. In order to understand the changes in ion release rate with composition and hence optimise Ag–PBG for specific applications, a detailed understanding of the atomic structure is required. First-principles molecular dynamics simulations using the density functional theory (DFT-MD) offer a parameter-free methodology accounting for the electronic structure of the system. These techniques have been extensively used to model PBG7–9 and phosphosilicate glasses10–15 allowing for an atomiclevel structural characterization and calculation of electronic properties. Classical techniques based on empirical force fields16 allow for simulations of larger cell sizes, compared to the prohibitive computational expense of DFT-MD, thus facilitating a reliable statistical description of medium-range structure.17 Furthermore, classical MD allows for slower glass quench rates (B1–10 K ps1), which have been shown to improve the microscopic and macroscopic descriptions of reference amorphous systems.18,19 These approaches have been used to model PBG and phosphosilicate glasses and provided structures in good agreement with experiment.20–24 The use of both classical and first-principles methodologies gives us greater confidence in our results by immediately highlighting any dependence of the results on simulation methodology. In this work we characterise, via first-principles and classical techniques, the short- and medium-range structure of experimentally synthesised3,5 glass compositions with 50 mol% P2O5,

Phys. Chem. Chem. Phys., 2014, 16, 21135--21143 | 21135

View Article Online

Paper

PCCP

0–20 mol% Ag2O, and varying amounts of Na2O and CaO. This represents the first step to understanding dissolution phenomena in Ag–PBG, from a theoretical perspective, via the unravelling of compositional effects on bulk structure.

Published on 14 July 2014. Downloaded by University of Connecticut on 28/10/2014 18:30:36.

II. Computational methods Born–Oppenheimer molecular dynamics (BOMD) simulations were performed with the QUICKSTEP module25,26 in the CP2K code (development version 2.3.16).27 A dual basis set was used, in which the Kohn–Sham orbitals are expanded in an atomcentered Gaussian basis set while the electronic charge density was described using an auxiliary plane-wave basis set.28 Core electrons were described with the pseudopotential of Goedecker, Teter and Hutter (GTH)29–31 incorporating scalar-relativistic core corrections. Valence electrons were treated with the Perdew, Burke and Ernzerhof (PBE)32 gradient-corrected exchange–correlation functional in the double-z valence polarized (DZVP) basis set. The orbital transformation method33 was employed for an efficient wavefunction optimization. The compositional nomenclature used henceforth is P = P2O5, C = CaO, N = Na2O and A = Ag2O, followed by the percentage molar composition for each component. The internal energy of a structurally randomised P50C30N10A10 composition is converged to 3.81  105 Ry per atom at a plane-wave kinetic-energy cutoff of 700 Ry for charge density, which was therefore used in the calculations. For every self-consistent field loop, the electronic gradient was converged to 1  105 Hartree. The timestep for each dynamics step was set to 1 fs with the positions, velocities and Mulliken charges34 of all the atoms recorded for each step. The glass generation for the first-principles (f.p.) calculations followed a standard melt-quench protocol in which atoms were initially quasi-randomly inserted into cubic simulation boxes (periodic in three dimensions) to reproduce the experimental densities given in Table 1. Two compositions (P50C40N10 and P50C30N10A10) were studied. For P50C30N10A10, a linear extrapolation of the experimental densities for 107Ag and 109Ag (ref. 35) was used, to reflect 50 : 50 isomeric abundancies. The configurations were then evolved for 3 ps in the microcanonical (NVE) ensemble. Subsequently the systems were heated to 2500 K in the canonical (NVT) ensemble and equilibrated for 30 ps. Equilibration was checked by examining the ions’ actual and meansquare displacements. The systems were then cooled to 300 K at a rate of B24 K ps1 in the NVT ensemble. The protocol followed 15 ps NVT trajectories run consecutively at 2200 K, 1900 K, 1600 K,

Table 1 Glass compositions simulated. First-principles and classical methodologies labelled as (f.p.) and (c.) respectively

Glass code

x

Atoms

Cell (Å)

1300 K, 1000 K and 700 K. At 300 K the trajectory was evolved for 40 ps during which time structural and electronic data were sampled. The f.p. results in this paper are averaged over this time using every timestep. By checking the ions’ mean-square displacements at 300 K, we can confirm that the model is a solid. The classical (c.) simulations followed a similar protocol, at a reduced quench rate of B2 K ps1, using the formalcharge, polarizable force field shown in Table 2 to describe inter-atomic forces. This has previously been shown to accurately reproduce the structure of glasses in the system P2O5– CaO–Na2O8,23,24 and now we include an Ag–Os two-body potential taken from Woodley et al.37 Molecular dynamics were run using the DL_POLY code38 (version 2.20) using an Evans39 thermostat in the NVT ensemble. In addition to the two compositions studied using first-principles methods, two compositions with higher silver content (identical to those prepared experimentally by Valappil et al.3) were simulated in order to compare to all available experimental data on dissolution rate. Three independent 3000-atom models were prepared of each of the high-Ag compositions and the results presented are averaged over all three models. The dynamics timestep was set at 0.2 fs for an accurate description of the Os motion. The classical systems were heated to 2400 K and melted for 50 ps, after which the temperature was reduced in steps of 100 K, running each trajectory for 50 ps. Core–shell frictional damping was parameterized as 20 r c2 r 25 (where c is the core–shell damping coefficient which is linear in velocity) for all temperatures in the range 2400–1500 K, c2 = 15 at 1400 K, c2 = 10 at 1300 K and c2 = 5 for all temperatures in the range 1200–300 K. A density correction factor of 0.950 was applied to the trajectories run at 2400 K, 2300 K and 2200 K. Likewise, a factor of 0.975 was applied at 2100 K, 2000 K and 1900 K in order to mimic thermal expansion. Since C(Os–Os) is non-zero, the simple power law expression for the dispersive energy means that V(Os–Os) - N for low r. To avoid this unphysical possibility (there is a finite

Table 2 Formal charge shell-model force field used in this work, including Buckingham-type two-body, harmonic three-body and core–shell (cs) potentials r r

Vij ¼ Ae

 Cr6

i–j

A (eV)

r (Å)

C (eV Å6)

P–Os23 Os–Os41 Na–Os42 Ca–Os42 Ag–Os37

1020.0000 22764.30 56465.3453 2152.3566 962.197

0.343220 0.149000 0.193931 0.309227 0.300000

0.0300 27.88 0.00 0.099440 0.00

Vijk = 12k3b(y  y0)2

Density (g cm3)

i–j–k

k3b (eV rad2)

y0 (1)

Os–P–Os23 P–Os–P23

3.3588 7.6346

109.470000 141.179333

P50C40N10 (f.p.) P50C30N10A10 (f.p.)

0.00 0.10

294 298

15.9862 16.2615

2.59036 2.84635

P50C40N10 (c.) P50C30N10A10 (c.) P50C30N5A15 (c.) P50C30A20 (c.)

0.00 0.10 0.15 0.20

6008 6007 2999 2999

43.6907 44.3763 35.0487 34.9521

2.59036 2.84635 3.094 3.331

21136 | Phys. Chem. Chem. Phys., 2014, 16, 21135--21143

Vij = 12kcsr 2 i–j

kcs (eV Å2)

Oc (e)

Os (e)

Oc–Os41

74.92

+0.8482

2.8482

This journal is © the Owner Societies 2014

View Article Online

PCCP

Paper

Published on 14 July 2014. Downloaded by University of Connecticut on 28/10/2014 18:30:36.

probability that the system will gain sufficient kinetic energy to overcome the repulsive barrier) during the high temperature trajectories 2400–1500 K, the potential was substituted for one of A the form VðOs Os Þ ¼ B (ref. 40 and 42), below a cut-off distance r of r(Os–Os) = 0.952 Å. As for f.p., all classical data are averaged from the 300 K trajectory with output printed every 1 fs (to match f.p. despite the smaller timestep of 0.2 fs).

III. Results and discussion A.

Short-range order: network formers

The glass network for each composition is made up of interlinked PO4 tetrahedra. Fig. 1 displays the partial pair-correlation functions g(P–O)(r) between 1.2–2.0 Å (within the first coordination shell). Both methodologies distinguish between the longer phosphorus to ‘‘bridging’’ oxygen bonds (P–BO), with g(P–O)(r) peaks centered between 1.62–1.63 Å and the shorter phosphorus to ‘‘non-bridging’’ oxygen bonds (P–NBO), centered between 1.47–1.50 Å. The slight overestimation from DFT-MD results are likely due to GGA underbinding and hence overestimating bond lengths. Pickup et al.36 used neutron diffraction (ND) and 31 P magic angle spinning nuclear magnetic resonance (MAS NMR) spectroscopy to probe the structure of P50C40N10. The structural parameters were obtained via simulation of the reciprocal-space

Fig. 1

data and conversion of the results to real space via Fourier transform (to make comparison to the correlation function from ND). Their results give r(P–BO) = 1.60  0.02 Å and r(P–NBO) = 1.49  0.02 Å. The classical simulations reproduce these bond lengths well (see Table 3), as does the glass simulated by first-principles techniques. For the Ag–PBG composition P50C30N10A10 comparison is made to data from ND experiments,35 using samples enriched with both 107Ag and 109Ag. Experimentally, it was shown that for both silver isotopes, r(P–BO) = 1.60  0.01 Å and r(P–NBO) = 1.48  0.01 Å. From the current work, r(P–BO) is slightly overestimated by both methodologies, with an accurate representation of r(P–NBO) in P50C30N10A10 (c.). It can be seen from these results that the addition of Ag to the glass, even at high-Ag content, has no effect on r(P–BO) and r(P–NBO). The disorder parameters sij, measuring static and thermal disorder, from the pair functions for P–BO and P–NBO are a measure of the width of the relevant pair-correlation peak. In experimental work, sij are derived from Q-space simulation using eqn (1): pðQÞij ¼

  Nij wij sin QRij Q2 sij2 ; exp cj QRij 2

(1)

where p(Q)ij is the pair function in reciprocal space, Nij, Rij and sij are the coordination number, atomic separation and disorder parameter (referred to as sND ij from this point onwards), respectively, of atom i with respect to j. cj is the concentration of atom j and wij is the weighting factor, given by wij = 2cicjbibj if i a j and wij = ci2bi2 if i = j (where b represents the coherent scattering length). From the current work the full-width-halfmaximum (FWHM) values of the decomposed partial-pair correlation functions g(P–O) (f.p.) have been calculated. Making the assumption that g(P–NBO)(r) and g(P–BO)(r) are normally distributed (for the first coordination sphere) gives, pffiffiffiffiffiffiffiffiffiffiffi f:p: f.p. FWHMf:p: ij ¼ 2 2 ln 2sij , from which sij are derived. Comparison of the compositional trends in simulated sND (P–BO) f.p. and sND (P–NBO) from ND experimental work, with trends in sij values of g(P–BO) and g(P–NBO) from this work, is presented in Table 4. Experimentally, it is noted that P–BO disorder increases upon ND 10 mol% Ag-doping from sND (P–BO) = 0.048 Å to s(P–BO) = 0.060 Å. f.p. In our calculations, s(P–BO) increases from 0.048 Å in P50C40N10 (f.p.) to 0.052 Å in P50C30N10A10 (f.p.) in good agreement with ND. sND (P–NBO) increases from 0.036 Å to 0.060 Å

P–O partial pair-correlation functions.

Table 3 Selected partial pair-correlation peak distances (r(X–X) (Å)). Experimental data obtained via neutron (ND) and X-ray (XRD) diffraction, are given for comparison

P50C40N10

r(P–NBO) r(P–BO) r(O–O) r(P–P) r(Ca–O) r(Na–O) r(Ag–O)a a

P50C30N10A10 36

f.p.

c.

ND

1.50 1.63 2.56 2.94 2.35 2.38 —

1.48 1.62 2.54 3.02 2.33 2.36 —

1.49 1.60 2.52 2.93 2.34 2.33 —

     

0.02 0.02 0.02 0.02 0.02 0.02

35

f.p.

c.

ND & XRD

1.50 1.63 2.57 2.96 2.34 2.37 2.37

1.48 1.62 2.55 3.03 2.31 2.34 2.29

1.48  0.01 1.60  0.01 2.51  0.01 2.93  0.01 2.38  0.01 2.33  0.01 2.28, 2.51, 2.73  0.03

P50C30N5A15

P50C30A20

c.

c.

1.48 1.63 2.55 3.03 2.32 2.36 2.27

1.47 1.63 2.55 3.03 2.32 — 2.29

Further analysis of Ag–O average distances (%r(Ag–O)) given in Table 6.

This journal is © the Owner Societies 2014

Phys. Chem. Chem. Phys., 2014, 16, 21135--21143 | 21137

View Article Online

Paper

PCCP

Published on 14 July 2014. Downloaded by University of Connecticut on 28/10/2014 18:30:36.

Table 4 Full-width-half-maximum (FWHM) values from the decomposed partial-pair correlation functions g(P–O) (f.p.) and disorder parameters sf.p. (standard deviation under the assumption that g(P–BO) and g(P–NBO) ij are normally distributed). Disorder parameters sND from Q-space simulaij tion (see eqn (1)) of experimental ND diffraction data35,36 f.p. ND FWHMf.p. ij (Å) sij (Å) sij (Å)

i–j

Model

P–BO P–BO

P50C40N10 (f.p.) 0.114 P50C30N10A10 (f.p.) 0.123

0.048 0.052

0.048  0.0136 0.060  0.0135

P–NBO P50C40N10 (f.p.) 0.075 P–NBO P50C30N10A10 (f.p.) 0.076

0.032 0.032

0.036  0.0136 0.060  0.0135

with the inclusion of silver. No significance can be attributed to the 0.001 Å increase upon doping, for the associated g(P–NBO) FWHM values from this work. However, ND and sf.p. (P–NBO) are in reasonable agreement for P50C40N10. Increases in disorder parameters are linked to changes in local bonding environments and, in turn, linked to a disproportionation in the mediumrange structure35 (see sub-section 3.4). Fig. 2 displays the partial pair-correlation functions for O–O and P–P. It can be seen that there are no significant compositional or methodological dependencies on the O–O distributions. However, the second coordination shell at B3.2 Å, is stronger for first-principles results in conjunction with stronger peaks in the O–P–O (f.p.) angular distribution function (ADF) (see Fig. 3), when compared to classical results. These differences are most likely a feature of the parameterization of the force field used in this work. The first peak in g(P–P)(r) shifts from a peak-centered position of 2.94 Å in P50C40N10 (f.p.) to 2.96 Å in P50C30N10A10 (f.p.) and sharpens considerably. This is accompanied by an increase of approximately 4% in the peakcentered positions of P(P–O–P)(y) from 1251 in P50C40N10 (f.p.) to 1291 in P50C30N10A10 (f.p.), thus indicating a possibly significant structural change upon doping the glass via substitution of CaO for Ag2O. This is further discussed in Section 3.3, in relation to the coordination environment of Ag. The classical results do not corroborate these features due to the inclusion of a P–Os–P three-body potential, as previously discussed.23

Fig. 2

O–O and P–P partial pair-correlation functions.

21138 | Phys. Chem. Chem. Phys., 2014, 16, 21135--21143

Fig. 3

B.

O–P–O and P–O–P angular distribution functions.

Short-range order: network modifiers

The partial pair-correlation functions for the network modifiers (Me) Ca, Na and Ag, with respect to oxygen, are given in Fig. 4. The classical and first-principles distributions are in good agreement. For both methodologies, the normalised peak intensities for Ca–O and Na–O increase upon substitution of CaO for Ag2O, with the consequence of decreasing the FWHM. The ill-defined first minimum for g(Ag–O)(r) from P50C30N10A10 (f.p.) (2.7–3.2 Å), suggests a distorted first coordination shell (however the 298 atom model also leads to poorer statistical sampling and increased signal noise). This in turn narrows the bond length distributions for the other modifiers Ca and Na due to geometrical constraints. The peaks of the distributions g(Me–O)(r) from classical simulations are centered on bond lengths 0.05–0.15 Å below the peaks of respective compositions from first-principles results, consistent with comparison of a structure derived from a wellparameterized force field and an underbinding GGA calculation. As with the short-range structure of the network-forming species, the network-modifier bond lengths are in good agreement with available experimental data (see Table 3).

Fig. 4 Ca–O, Na–O and Ag–O partial pair-correlation functions.

This journal is © the Owner Societies 2014

View Article Online

Published on 14 July 2014. Downloaded by University of Connecticut on 28/10/2014 18:30:36.

PCCP

Paper

Fig. 5 O–Ca–O, O–Na–O and O–Ag–O angular distribution functions. Me–O cut-offs taken as first minimum of respective radial distribution functions.

Fig. 5 displays the ADFs for the three-body systems O–Me–O. For the O–Ca–O system there is a strong localisation of bond angles at B801 and a significant contribution of angles up to 1501, consistent with a distorted octahedral coordination, as found in other Ag-free PBG.23 As also indicated, by a large normalised peak intensity in g(Ca–O)(r) and first minimum located at 3 Å, the O–Ca–O distribution for P50C30N10A10 (f.p.) displays increased intensity at B801 when compared to P50C40N10 (f.p.) (a dopant-dependent feature also seen in classical results). This indicates that the Ca ions in this model bond in a more symmetric octahedron which is likely linked to the presence of two further modifier species in the composition. A secondary set of peaks of lower intensity in the range 501–651 are attributed to intra-tetrahedral environments, where both oxygens are bonded to a common phosphorus atom.23,43 The ADFs for O–Na–O show significant compositional dependence for f.p. models. P50C40N10 (f.p.) has peaks at 571, 751 and 1201 compared to the stronger octahedral coordination of P50C30N10A10 (f.p.) with peaks at 521, 841 and B1601. Results from the classical simulations show no significant compositional dependence although there is increased intra-tetrahedral contribution. The O–Ag–O distribution for P50C30N10A10 (f.p.) has three dominant peaks centered at 541, 841, 1181 with a shoulder feature at B1401, suggesting a distorted bonding environment (representative Ag environment shown in Fig. 6). The classical results show very good agreement, with a slightly smoother distribution of O–Ag–O bond-angles (likely a statistical effect coupled with a differing description of inter-atomic forces). Peaks from the latter are centered on 541 and 841 with a shoulder at B1201 (see Fig. 5). C.

Modifier coordination

Both compositions studied contain a variety of local environments. f.p. results commonly show more highly resolved angular contributions, which may be a product of poor statistical sampling, or may reflect the highly accurate description of

This journal is © the Owner Societies 2014

Fig. 6 Example Ag coordination environment in P50C30N10A10 (f.p.) (upper panel) and crystalline Ag2SO4 – space group Fddd, No. 7045 (lower panel). Common distorted octahedral geometry and in-plane intratetrahedral coordination.

inter-atomic forces. As previously shown from ADFs, Ca displays a distorted octahedral coordination environment and this is confirmed by coordination numbers (Z(Ca–O)) of 6.91 and 6.57 for P50C40N10 (f.p.) and P50C30N10A10 (f.p.) respectively (see Table 5). This 0.34 decrease in Z(Ca–O) further indicates increased octahedral symmetry upon 10 mol% Ag-doping. f.p. results show Z(Na–O) = 6.07 reducing to Z(Na–O) = 5.85 when the glass is doped. These results are similar to previous MD simulations of ultra-phosphate PBG compositions.23 For f.p.

Phys. Chem. Chem. Phys., 2014, 16, 21135--21143 | 21139

View Article Online

Paper

PCCP

Published on 14 July 2014. Downloaded by University of Connecticut on 28/10/2014 18:30:36.

Table 5 Modifier & dopant coordination environments for f.p. and c. simulations. Coordination numbers (Z) with Ca–O, Na–O and Ag–O cutoffs set at 3.22 (f.p.)/3.25 (c.) Å, 3.15 (f.p.)/3.15 (c.) Å and 3.20 (f.p.)/3.15 (c.) Å respectively (bond lengths at first minima of respective partial paircorrelation functions). Z decomposed into BO and NBO contributions

P50C40N10 (f.p.)

P50C30N10A10 (f.p.)

Atomic pair

Z

Z

Ca–O Ca–BO Ca–NBO

6.91 0.52 6.39

6.57 0.48 6.09

Na–O Na–BO Na–NBO

6.07 0.99 5.08

5.85 0.85 5.00

Ag–O Ag–BO Ag–NBO

— — —

5.42 0.82 4.60

P50C40N10 P50C30N10A10 P50C30N5A15 P50C30A20 (c.) (c.) (c.) (c.) Atomic pair Z

Z

Z

Z

Ca–O Ca–BO Ca–NBO

6.79 0.65 6.14

6.44 0.49 5.95

6.57 0.47 6.10

6.51 0.73 5.78

Na–O Na–BO Na–NBO

6.48 1.43 5.05

6.02 1.18 4.84

6.26 1.29 4.98

— — —

Ag–O Ag–BO Ag–NBO

— — —

5.54 0.98 4.56

5.62 1.04 4.58

5.71 1.05 4.65

Table 6 Silver to oxygen coordination numbers (Z) as a function of Ag–O cut-off for P50C30N10A10 (f.p.)

2.80 Å 2.90 Å 3.00 Å 3.10 Å 3.20 Å 3.30 Å 3.40 Å Atomic pair Z

Z

Z

Z

Z

Z

Z

Ag–O Ag–BO Ag–NBO

4.07 0.33 3.74

4.47 0.47 4.00

4.92 0.63 4.29

5.42 0.82 4.60

5.98 1.05 4.93

6.55 1.29 5.26

3.71 0.22 3.49

results, ZAg–O = 5.42 at a cut-off of 3.20 Å. As previously stated, there is an ill-defined minimum in g(Ag–O)(r) and the dependence of Z(Ag–O) on Ag–O cut-off (Ag–O(cut)) is given in Table 6, ranging from Z(Ag–O) = 3.71 at Ag–O(cut) = 2.80 Å to Z(Ag–O) = 6.55 at Ag–O(cut) = 3.40 Å. Several X-ray absorption Near Edge Structure (XANES) studies of 0–20 mol% Ag PBG, have shown the shape and position of the Ag K-edge spectra to be identical to the reference material Ag2SO4.3,35,44 From these studies, the authors have concluded that the Ag ion in Ag–PBG has a structural environment very similar to that of Ag2SO4. Fig. 6 shows an example Ag ion coordination environment from the P50C30N10A10 (f.p.) model along with the local coordination environment of AgI in Ag2SO4 in the upper and lower panels respectively.45 Both structures show a distorted octahedral symmetry, providing further evidence of the similarities in Ag local ordering.

21140 | Phys. Chem. Chem. Phys., 2014, 16, 21135--21143

Moss et al.35 have previously conducted neutron diffraction with isotopic substitution (NDIS) using P50C30N10A10 samples enriched with 107Ag and 109Ag isotopes. This allows for the simplification of complex correlations46 (assuming each sample has identical composition and structure) by the application of difference function methods.47,48 They showed35 that the Ag–O correlation has three components in the first coordination shell at 2.28 Å, 2.51 Å and 2.73 Å with coordination numbers of 2.1, 2.7 and 1.1 respectively. g(Ag–O)(r) from the current work (see Table 4 for peak position) displays a unimodal distribution within the first coordination sphere, for both methodologies, with no evidence for multiple resolved peaks. Table 7 shows a direct comparison of the Ag–O mean c. bond lengths for both methodologies, r%f.p. (Ag–O) and r%(Ag–O), in P50C30N10A10, as a function of Ag–O(cut) for the first coordina35 tion sphere. The experimental value of r%NDIS (Ag–O) = 2.47  0.07 Å is derived using the published bond lengths weighted accordc. ing to coordination number. r%f.p. (Ag–O) and r%(Ag–O) provide a more direct comparison to experiment, compared to the peak centered position of g(Ag–O)(r), due to the fact that the partial paircorrelation does not decay to zero at the limit of the first coordination sphere, resulting in the deviation of peak centered values from the mean bond length. r%f.p. %c. (Ag–O) and r (Ag–O) range from 2.46–2.61 Å and 2.41–2.64 Å respectively, across the range Ag–O(cut) = 2.90–3.30 Å, in good agreement with experiment. Fig. 6 facilitates a more detailed assignment of peaks in the P50C30N10A10 (f.p.) O–Ag–O ADF from Fig. 5. The peak at 541 corresponds to intra-tetrahedral bonding as evidenced in the example environment shown. This distorted octahedral geometry also shows further in-plane angles at 771, 1351 and 1471 contributing to the peak at 841 and shoulder feature at B1401 in the ADF. Furthermore, Ag–O first-coordination environments with 5-coordinated distorted trigonal bipyramidal geometries were observed, constituting contributions to the ADF peaks centered on 841 and 1181. Z(Me–O) values from classical simulations (c.) of P50C30N10A10 are given in the lower half of Table 5 with amended cut-offs to reflect differing first minima in g(Me–O)(r) (c.). Direct comparison to f.p. results shows good agreement for Z(Ca–O) and Z(Na–O) with similar decompositions for BO and NBO. Upon 10 mol% Ag-doping Z(Ca–O) (c.) and Z(Na–O) (c.) decrease by 0.35 and 0.46 respectively, compared to decreases of 0.34 (f.p.) and 0.22 (f.p.) respectively. Z(Ag–O) = 5.54 (c.) at Ag–O(cut) = 3.15 Å, which compared to similar cut-off values for (f.p.) results given in Table 6, and 5.62 and 5.71 for high-Ag content compositions, is in good agreement. Table 7 Mean silver to oxygen bond lengths (%r(Ag–O)) in first coordination sphere as a function of Ag–O cut-off for P50C30N10A10 (f.p.) and P50C30N10A10 (c.). Average bond length from experiment r%NDIS (Ag–O) = 2.47  0.07 Å35

2.90 Å 2.95 Å 3.00 Å 3.05 Å 3.10 Å 3.15 Å 3.20 Å 3.25 Å 3.30 Å r%f.p. (Ag–O)

2.46

2.47

2.49

2.50

2.53

2.54

2.57

2.58

2.61

(Å) r%(c. Ag–O) (Å)

2.41

2.44

2.47

2.50

2.53

2.55

2.58

2.61

2.64

This journal is © the Owner Societies 2014

View Article Online

PCCP

Published on 14 July 2014. Downloaded by University of Connecticut on 28/10/2014 18:30:36.

D.

Paper

Medium-range order

The larger simulation models used in classical methods facilitate a more statistically sound analysis of the medium-range structure, providing insight into compositionally dependent structural trends at these length scales. Thus, only results from the classical simulations are presented in this section. Moss et al.35 have inferred that a disproportionation in the mediumrange Qn distribution occurs upon CaO substitution with Ag2O for the compositions studied here. The authors, using a simple bond-order model, show that an increase in s(P–BO) and static disorder of P–NBO implies a disproportionation in the network structure. Valappil et al.3 explicitly showed, via 31P MAS NMR spectral analysis, that the substitution of Na2O with Ag2O, for the compositions (P2O5)0.50(CaO)0.30(Na2O)0.20–x(Ag2O)x (x = 0.00, 0.10, 0.15 and 0.20), leads to a disproportionation of Q2 units approximating the relation 2Q2 - Q1 + Q3, for x Z 0.10. The phosphorus Qn distributions from classical simulations are given in Table 8. P50C40N10 (c.) comprises 79.7% Q2, 9.9% Q1 and 9.7% Q3 from our simulations compared to 96% Q2 and 4% Q1 from experiment.36 Similar levels of agreement with experiment are found for classical MD derived phosphate,23 phosphosilicate20 and sodium silicate42,49 glasses. The potential reasons for discrepancies between theoretical and experimental Qn distributions are wide ranging and include modal assumptions when deconvoluting experimental NMR signals. From the theory side, Stebbins50 found that both an increase in glass transition temperature (known to be higher in theoretical simulations than experiment) and higher field strength modifiers caused the disproportionation reactions, 2Qn - Qn1 + Qn+1 (for n = 1, 2 and 3), in silicate glasses. The discrepancies in Q1, Q2 and Q3 proportions for P50C40N10 (c.) from the current work (compared to experiment) may relate to fictive temperatures that are in excess of experiment (further limitations with MD methods are discussed by Pota et al.51). In any case, we are interested in the trends of the Qn distribution on Ag-doping, and the absolute values are less important. For P50C30N10A10 (c.), Q2 reduces to 78% with concomitant rises in Q1 and Q3, approximating the relation found by Valappil and co-workers,3 but we do not observe this disproportionation occurring at 410 mol% Ag2O (Table 8). Network connectivity (NC), defined as the mean number of BO within the first coordination shell of phosphorus summed over all P atoms, is 1.99–2.00 for all simulated compositions. This is very close to the value of 2.0 expected theoretically and found experimentally.36

Table 8 Qn species distribution (%) (for phosphorus with respect to oxygen) and network connectivity (NC). Experimental data (expt.) from fitting to 31P MAS NMR spectra36 0

P50C40N10 (expt.) P50C40N10 (c.) P50C30N10A10 (c.) P50C30N5A15 (c.) P50C20A20 (c.)

1

2

3

Q

Q

Q

Q

0.0 0.6 0.4 0.2 0.3

4.0 9.9 10.9 9.0 9.3

96.0 79.7 78.0 82.0 81.1

0.0 9.7 10.6 8.8 9.3

This journal is © the Owner Societies 2014

Q

4

0.0 0.1 0.1 0.0 0.1

NC 1.96 1.99 1.99 1.99 2.00

Table 9 Mean Mulliken charges for all atomic species X (Q(X)) in P50C40N10 (f.p.) and P50C30N10A10 (f.p.)

P50C40N10 (f.p.)

P50C30N10A10 (f.p.)

Q(P) Q(O)

+0.760 0.459

+0.766 0.453

Q(BO) Q(NBO)

0.384 0.497

0.384 0.488

Q(Ca) Q(Na) Q(Ag)

+1.146 +0.810 —

+1.130 +0.791 +0.540

E.

Atomic charges

The mean atomic charges obtained by Mulliken population analysis are given in Table 9. These are summed over all timesteps of the 300 K trajectory from f.p. models. No significant changes are noted upon Ag-doping with net phosphorus charges of +0.760 and +0.766 for P50C40N10 (f.p.) and P50C30N10A10 (f.p.) respectively. It should be noted that Mulliken charges are calculated by determining the electron population of each atom as defined by the basis functions and absolute values are therefore dependent on the choice of basis set. Hence, although the partition of the electronic density is somewhat arbitrary, we use Mulliken charge analysis for direct comparison to previous work, and for its relative ease of implementation. We are not attempting any detailed quantitative analysis of the atomic charges. The electronic states and chemical bonding in phosphate glass have previously been studied by Kowada et al.52 using the DV-Xa cluster method. For the charge-balanced P4O10 cluster the net charge on phosphorus was found to be +1.31 using Mulliken population analysis. In the context of the other atoms in this study, the phosphorus atoms have a covalent interaction with the surrounding oxygens. The mean oxygen atomic charges have been deconvoluted into BO and NBO, with the NBO displaying increased ionic character with net charges of 0.497 and 0.488 for PBG and Ag–PBG respectively. Mean BO charges are found to be +0.113 and +0.104 higher due to the presence of two covalently bonded phosphorus atoms in their first coordination sphere. These changes are in good agreement with the relation Q(BO) = Q(NBO) + 0.12 for the cluster P4O10.52 All modifier ions show increased ionic character compared to the network formers. XANES analysis conducted by Ahmed et al.44 has shown silver to be in the AgI oxidation state for a series of Ag–PBG. The effective ionic radii for six-coordinate AgI and NaI are 1.15 and 1.02 Å respectively,53 suggestive of a higher effective charge for NaI compared to AgI, as corroborated by Q(Na) = +0.791 and Q(Ag) = +0.540 for P50C30N10A10 (f.p.).

IV. Discussion and conclusions We have investigated the structures of biomedically relevant glasses with 50 mol% P2O5, 0–20 mol% Ag2O and varying amounts of Na2O and CaO, using first-principles (f.p.) and classical (c.) molecular dynamics simulations. Good agreement is found between the results from c. simulations (using a force

Phys. Chem. Chem. Phys., 2014, 16, 21135--21143 | 21141

View Article Online

Published on 14 July 2014. Downloaded by University of Connecticut on 28/10/2014 18:30:36.

Paper

field that was previously developed within our group23) and the high level f.p. reference calculations. The c. methodology outperforms f.p. in reproducing some key structural features, when compared to experimental data, such as P–O bond distances. The Ag ion is noted to have little effect on shortrange order. Ag itself has a distorted octahedral and trigonal bipyramidal structure with an ill-defined first coordination shell where the Ag–O coordination number is approximately 5.5 and the bond length is B2.5 Å, implying that Ag is acting as a network modifier. At low Ag content, and for Ag2O 2 CaO substitution, there is experimental evidence for a disproportionation in the medium-range Qn distribution for phosphorus upon doping, following the relation 2Q2 - Q1 + Q3, which we also find in simulation. An increase in P–BO bond disorder (as seen in experiment35,36) from f.p. results on 0 and 10 mol% Ag2O compositions, further supports this feature. These structural changes are of likely consequence in the observed decrease in dissolution rates with increased Ag mol%. However, we do not see this disproportionation at higher Ag contents when Ag2O is substituted for Na2O. In fact, it is rather difficult to disentangle the precise effect of Ag on the glass dissolution because of the different substitutions employed in past experiments. We cannot say with certainty whether the composition dependence of the disproportionation reaction is due to differing Ag content, or to differences in the moiety for which Ag2O is substituted. What is perfectly clear from experiment is that when Ag2O is substituted for Na2O, the dissolution rate goes down. Our simulations show clearly that this is due to differences in the local bonding. Ag has a field strength intermediate between Na and Ca, which is reflected in the intermediate number of bridging oxygen atoms to which Ag is bonded (Table 5). The experimental implication is that Ag is bonded to more parts of the phosphate network than Na, and that its substitution for Na strengthens the glass network through greater cross-linking. Although we show that Ag has a lower coordination number than Na, it seems that this is less important for the dissolution than the value of its field strength. Although it is doubtful that 50 mol% P2O5 glasses have short phosphate chains in the way that lower-phosphate content glasses do, the effect is similar to Ca bonding together more phosphate chains than Na in lowphosphate glasses, as we have previously identified.24

Acknowledgements The authors would like to thank the Engineering and Physical Sciences Research Council (United Kingdom) EPSRC(GB) (EP/ J008095) for funding. Via our membership of the UK’s HPC Materials Chemistry Consortium, which is funded by EPSRC (EP/F067496), this work made use of the facilities of HECToR, the UK’s national high-performance computing service, which is provided by UoE HPCx Ltd at the University of Edinburgh, Cray Inc. and NAG Ltd, and funded by the Office of Science and

21142 | Phys. Chem. Chem. Phys., 2014, 16, 21135--21143

PCCP

Technology through EPSRC’s High End Computing Programme. We are grateful to Devis Di Tommaso for helpful discussions.

References 1 R. Sheridan, P. J. Doherty, T. Gilchrist and D. Healy, J. Mater. Sci.: Mater. Med., 1995, 6, 853–856. 2 T. Gilchrist, D. M. Healy and C. Drake, Biomaterials, 1991, 12, 76–78. 3 S. P. Valappil, D. M. Pickup, D. L. Carroll, C. K. Hope, J. Pratten, R. J. Newport, M. E. Smith, M. Wilson and J. C. Knowles, Antimicrob. Agents Chemother., 2007, 51, 4453–4461. 4 A. M. Mulligan, M. Wilson and J. C. Knowles, Biomaterials, 2003, 24, 1797–1807. 5 I. Ahmed, D. Ready, M. Wilson and J. C. Knowles, J. Biomed. Mater. Res., Part A, 2006, 79A, 618–626. 6 C. P. Randall, L. B. Oyama, J. M. Bostock, I. Chopra and A. J. Oneill, J. Antimicrob. Chemother., 2013, 68, 131–138. 7 E. Tang, D. Di Tommaso and N. H. de Leeuw, Adv. Eng. Mater., 2010, 12, B331–B338. 8 D. Di Tommaso, R. I. Ainsworth, E. Tang and N. H. de Leeuw, J. Mater. Chem. B, 2013, 1, 5054. 9 J. K. Christie, R. I. Ainsworth and N. H. de Leeuw, Biomaterials, 2014, 35, 6164–6171. 10 A. Tilocca, N. H. de Leeuw and A. N. Cormack, Phys. Rev. B: Condens. Matter Mater. Phys., 2006, 73, 104209. 11 A. Tilocca, Phys. Rev. B: Condens. Matter Mater. Phys., 2007, 76, 224202. 12 A. Tilocca, A. N. Cormack and N. H. de Leeuw, Faraday Discuss., 2007, 136, 45–55. 13 J. K. Christie and A. Tilocca, Adv. Eng. Mater., 2010, 12, B326–B330. 14 A. Tilocca, J. Mater. Chem., 2010, 20, 6848–6858. 15 J. K. Christie, A. Pedone, M. C. Menziani and A. Tilocca, J. Phys. Chem. B, 2011, 115, 2038–2045. 16 A. Pedone, J. Phys. Chem. C, 2009, 113, 20773–20784. 17 A. Tilocca, J. Chem. Phys., 2008, 129, 084504. 18 K. Vollmayr, W. Kob and K. Binder, J. Chem. Phys., 1996, 105, 4714–4728. 19 K. Vollmayr, W. Kob and K. Binder, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 15808–15827. 20 A. Tilocca, A. N. Cormack and N. H. de Leeuw, Chem. Mater., 2007, 19, 95–103. 21 A. Tilocca, Proc. R. Soc. London, Ser. A, 2009, 465, 1003–1027. 22 J. K. Christie and A. Tilocca, Chem. Mater., 2010, 22, 3725–3734. 23 R. I. Ainsworth, D. Di Tommaso, J. K. Christie and N. H. de Leeuw, J. Chem. Phys., 2012, 137, 234502. 24 J. K. Christie, R. I. Ainsworth, D. Di Tommaso and N. H. de Leeuw, J. Phys. Chem. B, 2013, 117, 10652–10657. 25 M. Krack, M. Parrinello, in High Performance Computing in Chemistry, ed. J. Grotendorst, 2004, vol. 25, pp. 29–51. 26 J. VandeVondele, M. Krack, F. Mohamed, M. Parrinello, T. Chassaing and J. Hutter, Comput. Phys. Commun., 2005, 167, 103–128.

This journal is © the Owner Societies 2014

View Article Online

Published on 14 July 2014. Downloaded by University of Connecticut on 28/10/2014 18:30:36.

PCCP

27 CP2K developers group: http://www.cp2k.org; 2000–2013. 28 G. Lippert, J. Hutter and M. Parrinello, Mol. Phys., 1997, 92, 477–488. 29 S. Goedecker, M. Teter and J. Hutter, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 1703–1710. 30 C. Hartwigsen, S. Goedecker and J. Hutter, Phys. Rev. B: Condens. Matter Mater. Phys., 1998, 58, 3641–3662. 31 M. Krack, Theor. Chem. Acc., 2005, 114, 145–152. 32 J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868. 33 J. VandeVondele and J. Hutter, J. Chem. Phys., 2003, 118, 4365. 34 R. S. Mulliken, J. Chem. Phys., 1955, 23, 1833–1840. 35 R. M. Moss, D. M. Pickup, I. Ahmed, J. C. Knowles, M. E. Smith and R. J. Newport, Adv. Funct. Mater., 2008, 18, 634–639. 36 D. M. Pickup, I. Ahmed, P. Guerry, J. C. Knowles, M. E. Smith and R. J. Newport, J. Phys.: Condens. Matter, 2007, 19, 415116. 37 S. M. Woodley, P. D. Battle, J. D. Gale and C. R. A. Catlow, Phys. Chem. Chem. Phys., 1999, 1, 2535–2542. 38 W. Smith and T. R. Forester, J. Mol. Graphics, 1996, 14, 136–141. 39 D. J. Evans and G. P. Morriss, Comput. Phys. Rep., 1984, 1, 297–343. 40 V. A. Bakaev and W. A. Steele, J. Chem. Phys., 1999, 111, 9803–9812.

This journal is © the Owner Societies 2014

Paper

41 M. J. Sanders, M. Leslie and C. R. A. Catlow, J. Chem. Soc., Chem. Commun., 1984, 19, 1271. 42 A. Tilocca, N. H. de Leeuw and A. N. Cormack, Phys. Rev. B: Condens. Matter Mater. Phys., 2006, 73, 104209. 43 J. K. Christie, J. Malik and A. Tilocca, Phys. Chem. Chem. Phys., 2011, 13, 17749. 44 I. Ahmed, E. A. Abou Neel, S. P. Valappil, S. N. Nazhat, D. M. Pickup, D. Carta, D. L. Carroll, R. J. Newport, M. E. Smith and J. C. Knowles, J. Mater. Sci., 2007, 42, 9827–9835. ¨pke and A. Illguth, 45 B. N. Mehrotra, T. Hahn, W. Eysel, H. Ro Neues Jahrb. Mineral., Monatsh., 1978, 408–421. 46 J. C. Wasse, I. Petri and P. S. Salmon, J. Phys.: Condens. Matter, 2001, 13, 6165–6176. 47 I. T. Penfold and P. S. Salmon, Phys. Rev. Lett., 1990, 64, 2164–2167. 48 P. S. Salmon, S. Q. Xin and H. E. Fischer, Phys. Rev. B: Condens. Matter Mater. Phys., 1998, 58, 6115–6123. 49 J. Du and A. N. Cormack, J. Non-Cryst. Solids, 2004, 349, 66–79. 50 J. F. Stebbins, J. Non-Cryst. Solids, 1988, 106, 359–369. 51 M. Pota, A. Pedone, G. Malavasi, C. Durante, M. Cocchi and M. C. Menziani, Comput. Mater. Sci., 2010, 47, 739–751. 52 Y. Kowada, H. Adachi and T. Minami, J. Phys. Chem., 1993, 97, 8989–8992. 53 R. D. Shannon, Acta Crystallogr., Sect. A: Cryst. Phys., Diffr., Theor. Gen. Cryst., 1976, 32, 751–767.

Phys. Chem. Chem. Phys., 2014, 16, 21135--21143 | 21143