A Computational Study of Magnesium Incorporation in the Bulk and

Article pubs.acs.org/Langmuir

A Computational Study of Magnesium Incorporation in the Bulk and Surfaces of Hydroxyapatite Neyvis Almora-Barrios,*,†,‡ Ricardo Grau-Crespo,† and Nora H. de Leeuw*,† †

Department of Chemistry, University College London, 20 Gordon Street, London WC1H 0AJ, United Kingdom ABSTRACT: We have used a combination of static lattice energy minimization and molecular dynamics simulations to investigate the thermodynamics of Mg incorporation into the bulk and hydrated surfaces of hydroxyapatite (HA). In agreement with recent experimental and theoretical work, our simulations show that the incorporation of low levels of Mg in the Ca (II) site is preferred with respect to incorporation in Ca (I) sites. However, we predict that Mg in the HA bulk material is metastable both with respect to the Mg/Ca exchange with aqueous solution and with respect to separation into bulk phases of magnesium phosphate [Mg3(PO4)2] and magnesium hydroxide [Mg(OH)2]. This finding suggests that Mg siting in the HA bulk is at least partially controlled by kinetics rather than by thermodynamics during crystal growth, which can explain the discrepancies found in the literature about the preferential substitution site. Finally, we found that Mg incorporation from solution into the hydrated surfaces, rather than the bulk material, is energetically favorable, in particular in the (011̅0) plane where cation sites are exposed to solution, thereby enabling the favorable interaction of Mg with water.

1. INTRODUCTION Defects and impurities in biological apatite, Ca10(PO4)6(F, Cl, OH)2, can be identified as substitutional defects or as discrete, extraneous crystalline phases. Ions that may be incorporated into the apatite structure, either intentionally or unintentionally, include carbonate ions1,2 (substituting hydroxy or phosphate groups), fluoride ions3 (substituting hydroxy groups), silicon (substituting phosphorus to form silicate groups),4 and magnesium and strontium ions (substituting calcium).5,6 In the case of hydroxyapatite (HA) forming bone, the incorporated impurities are released into, or taken up from, the surrounding body fluids during the normal bone remodelling processes. The content of Mg in HA is low (between ∼0.4 and 1.5 wt %),7 but it appears to play an important function in bone metabolism. Mg deficiency gives rise to decreased bone mass, owing to fewer osteoblasts and increased levels of osteoclasts.8 A recent paper has discussed the distribution of Mg substitution impurities over the two types of cation sites in HA, normally occupied by calcium [Ca(I) and Ca(II)], based on both theoretical and experimental evidence.9 There, we reported our prediction, based on interatomic potential calculations, that Mg preferentially substitutes at Ca(II) sites, which was found to be in agreement with NMR and XAS experiments and with a previous experimental study.10 However, other previous work had reported preferential Mg occupancy of the Ca(I) sites.11−13 The reasons for these discrepancies are not well-understood yet; they could be due either to lack of resolution in experimental measurements or to true variability of the distribution from sample to sample, depending on the synthesis route. The latter scenario would © 2013 American Chemical Society

imply that thermodynamics does not control the cation distribution in natural samples. In order to shed further light onto the fate of substitutional Mg defects in HA, in this paper we focus on the thermodynamics of Mg2+ incorporation in the HA bulk and surfaces, using atomistic simulation techniques in combination with statistical mechanics. In particular, we will discuss the stability of Mg in HA based on the equilibrium with aqueous solution of Ca and Mg ions and with respect to separation into competing Mg-bearing phases. We will also study the segregation of Mg from the hydrated (0001) and (011̅0) surfaces into the bulk material.

2. METHODOLOGY Natural hydroxyapatite has a hexagonal crystal structure with space group P63/m. The a and c parameters for hydroxyapatite are 9.425 Å and 6.884 Å, respectively.14 Four calcium ions in the unit cell are usually labeled as columnar calcium [or Ca(I)] because they form single atomic columns perpendicular to the basal plane. The other six out of ten calcium ions in the unit cell, labeled Ca(II), form hexagonal channels along the c direction of the structure, at the vertices of triangles around each hydroxy group. Two OH positions per unit cell are available in each channel, with a 0.5 occupancy of each position.15 In order to translate this structure into a model with full occupancies for the hydroxy groups, as is required to carry out the calculations, we have assigned alternate 0 and 1 occupancies to these sites along the hydroxy channels in the c direction, thus changing the space group of the hydroxyapatite unit cell from P63/m to P63,16 although this P63 cell is not yet suitable for realistic calculations either. Since there is only Received: February 1, 2013 Revised: April 12, 2013 Published: April 15, 2013 5851

dx.doi.org/10.1021/la400422d | Langmuir 2013, 29, 5851−5856

Langmuir

Article

Figure 1. HA structure with (P21) symmetry. Ca(II) atoms are at the vertices of triangles around each hydroxy group. (Ca = green, O = red, P = purple, H = white). one hydroxy channel per unit cell in the P63/m or P63 structures, all the hydroxy groups would be oriented in the same direction, creating a net electric polarization per unit cell. In the real HA structure, disorder in the relative orientation of the OH groups in the parallel channels means that electric polarization is not present in the material. Therefore, in order to create a more realistic structure, we have doubled the cell in the b direction and assigned opposite orientations to the OH groups in the two channels of the supercell. This antiparallel orientation of the hydroxy groups coincides with the arrangement of the OH groups in pure, synthetic HA, which crystallizes as an ordered monoclinic structure (P21) with a double unit cell compared to the hexagonal structure (Figure 1). In addition to the bulk material, we have also investigated two surfaces with low Miller indices, (0001) and (011̅0), which are modeled as periodic slabs of ∼20 Å thickness perpendicular to the surface. The slabs were separated from their images in the neighboring cells by a region of approximately 50 Å, which was filled with water molecules at a density of 0.99 g mL−1. Tests of the convergence of the total energy and properties with respect to the slab thickness showed that the surface energies are well-converged (within 0.01 J m−1). We have focused on the (0001) and (011̅0) surfaces because they are prominent in the morphology of the apatite mineral platelets in natural bone, dentin, and tooth enamel, and there is experimental evidence that these faces act as the binding sites for many ionic species, including small molecules, polymers, and anionically modified cell surfaces.17 The bulk lattice energies were calculated with the general utility lattice program (GULP),18−20 employing interatomic potential simulation techniques based on the Born model of solids,21 where it is assumed that the ions in the crystal interact via long-range electrostatic forces and short-range forces, including both the repulsions and van der Waals attractions between neighboring electron charge clouds. The interatomic potential model for hydroxyapatite is taken from previous work by de Leeuw22 and includes electronic polarizability via the shell model of Dick and Overhauser,23 where each polarizable ion, here oxygen, is represented by a core and a massless shell connected by a spring. The polarizability of the oxygen ions is then determined by the spring constant and the charges of the core and shell. Angle-dependent forces are included to allow directionality of bonding, for example in the phosphate groups. In order to make

geometry predictions, the lattice energies are minimized with respect to the structural parameters. All structures reported are the result of constant pressure energy minimizations, with an external pressure set to zero, where not only the ionic positions but also the cell parameters are allowed to vary to find the energy minimum. The distribution of Mg ions in the hydroxyapatite supercell was investigated using the site occupancy disorder (SOD) program,24,25 which generates the complete configurational space for each composition in a supercell of the structure before extracting the subspace of symmetrically independent configurations. The relative energies of the independent configurations are then calculated using GULP. The probability of occurrence of a given independent configuration m, with energy Em and degeneracy Ωm (the number of times it appears in the full configurational space), is given by Pm =

Ωm exp(− Em /RT ) Z

(1)

where Z is the partition sum, which guarantees that the sum of all the probabilities equals one; R is the gas constant, and T is the temperature. We take into account the dynamic nature of the surfaces in contact with water by using molecular dynamics simulations (MD), which were performed using the DL_POLY2 code.26 The potential model for water and for the interaction of water with the HA surfaces were taken from de Leeuw and co-workers.27−29 All surface simulations were performed at constant volume (cell parameters obtained from constant pressure calculations in the bulk) and temperature (310 K), using the Nosé−Hoover algorithm for the thermostat with a relaxation time of 0.5 ps.30,31 The integration of the equations of motion was performed using the Verlet leapfrog scheme.32 The total simulation time for each job was 0.5 ns, including 100 ps of equilibration.

3. RESULTS AND DISCUSSION 3.1. Magnesium Distribution in the Bulk. First, we review briefly our results of the Mg distribution over the bulk crystal sites of HA, which were originally presented in a recent paper,9 for comparison with experimental results and the surface calculations on ion exchange from solution in this work. 5852

dx.doi.org/10.1021/la400422d | Langmuir 2013, 29, 5851−5856

Langmuir

Article

We have considered the substitution of Mg in the different Ca positions in a 1 × 2 × 2 supercell with initial composition Ca40(PO4)24(OH)8. In this supercell model, the nearest Mg− Mg separation for a single substitution is 9.6 Å, along the c axis. Five inequivalent configurations were investigated for the isolated Mg atom in the different Ca positions, whereas 105 configurations were identified for two substituted Mg ions and 1235 configurations for three Mg substitutions, giving impurity concentrations of 2.5, 5, and 7.5 mol %, respectively. Of the five configurations for the 2.5% Mg composition (1 Mg per cell), two are of the Ca (I) type and three of the Ca(II) type. Inequivalent positions of the same type differ slightly from each other in their distances to the neighboring O atoms. The calculated energies indicated a preference for the Mg cation to be incorporated at a Ca(II) site, in agreement with both ab initio calculations and 43Ca solid state NMR analysis.9 In accordance with our calculations, the most stable configuration of Mg in Ca(I) was less energetically favorable by almost 34 kJ mol−1 compared to Mg in the Ca(II) position. The results of our calculations for 5% Mg still show a preference for magnesium cations to be incorporated in the Ca(II) sites. In the most stable configuration, the two Mg ions are located at alternating vertices of different adjacent triangles around the hydroxy groups in one channel. In the case of 7.5% Mg, (3 Mg per cell), the most stable configurations are again those where Mg ions are located in Ca(II) sites, while the mixed Mg(II) and Mg(I) configurations were higher in energy than the Mg(II)-only configurations (by ∼24 kJ mol−1). In summary, our results show that substitution of Mg for Ca in the bulk HA material should occur preferentially at the Ca(II) sites, at least for substitutions at low concentrations. 3.2. Stability of Mg in Hydroxyapatite with Respect to Phase Separation. We have now extended our previous work, which was a simple comparison of relative stabilities between two bulk sites without taking into account potential phase separation or alternative external environments (e.g., the surrounding aqueous solution). Here, we evaluate the stability of Mg in HA with respect to separation into other bulk phases. Because pure Mg hydroxyapatite is unknown, we calculate the enthalpy of mixing between pure Ca-hydroxyapatite plus a combination of magnesium phosphate [Mg3(PO4)2] and magnesium hydroxide [Mg(OH)2], as opposed to formation of the Mg-substituted hydroxyapatite. We have chosen to study the separation of Mg from HA, assuming a range of conditions under which HA itself is stable with respect to decomposition, as these are the conditions of relevance for applications as biomaterials. Therefore, although it is known experimentally that at high temperatures Mg should incorporate into the decomposition products of HA and in particular in tricalcium phosphate phases, these competing phases will not be considered in the present study

H (x ) =

(3)

In this expression, N is the number of formula units in the simulation cell. Figure 2 shows the enthalpies of mixing as a function of Mg concentration, calculated using eq 2. All values are positive and

Figure 2. Calculated enthalpies of mixing between hydroxyapatite and magnesium phosphate/magnesium hydroxide under ambient conditions.

exhibit an almost linear increase in the dilute limit. We have fitted our data to a quadratic curve in the form ΔHmix = Wx(1 − x)

(4)

where W is interpreted as the solution energy of Mg in HA from the pure phosphate and hydroxide. At T = 300 K, we obtained W = 82.3 kJ/mol, which is very high compared to RT = 2.5 kJ/mol. The maximum concentration of Mg that is stable with respect to separation into magnesium phosphate and magnesium hydroxide phases can be estimated from the minimum of the mixing free energy, which at very low impurity concentrations (and therefore ideal disorder) can easily be expressed as a function of the fraction x of substituted sites: Gmix (x) = Wx(1 − x) + RT[x ln x + (1 − x) ln(1 − x)] (5)

The free energy remains positive even after considering the negative configurational entropy contribution, except when x ≤ e−W/RT, which is 10−18 at T = 300 K. Any Mg concentration higher than that (and therefore, any significant Mg content) is only metastable with respect to separation into a magnesium phosphate and magnesium hydroxide phase. The immiscibility of Mg in hydroxyapatite is not unexpected if we consider the significant size difference between the two ions (r[Ca2+] = 1.00 Å and r[Mg2+] = 0.72 Å in six-coordinated environments).33 3.3. Stability of Mg in Hydroxyapatite with Respect to Exchange with Aqueous Solution. Phase separation is not the only limiting factor for Mg incorporation in the bulk, as the equilibrium concentration of the impurity will also depend on the relative concentrations of the ions in aqueous solution. We have calculated the average energies of isolated Mg2+ and Ca2+ cations in water, also using molecular dynamics at 300 K. The simulations give very reasonable results both in terms of energy and of the structure of the water layer around the ions. The 2+ calculated energy difference E[Mg2+ (aq)] − E[Ca(aq)] is equivalent to the difference in hydration energies between the two cations, 2+ ΔHhyd = [Mg2+ (aq)] − Hhyd[Ca(aq)], which is found experimentally

ΔHmix = H(x) − (1 − x)H(0) − xH {3[Mg3(PO4 ]2 ) + Mg(OH)2 }

1 ∑m EmΩm exp( −Em /RT ) N ∑m Ωm exp( −Em /RT )

(2)

where H(0) is the calculated energy per formula unit of pure Ca−HA and H(x) is the enthalpy per formula unit of the solid solution with a fraction x of Ca sites substituted by Mg, which is obtained from a Boltzmann-weighted average of the configuration energies Em (ignoring pressure and vibrational contributions to the enthalpy) for that composition: 5853

dx.doi.org/10.1021/la400422d | Langmuir 2013, 29, 5851−5856

Langmuir

Article

to be −344 kJ mol−1 at ambient temperature.34 Our theoretical result at 300 K is ΔHhyd = −344.9 kJ mol−1, in excellent agreement with the experiment. The very negative value indicates that the interaction with water is much more energetically favorable for Mg2+ than for Ca2+ cations. We can now consider the equilibrium between the solid phase and the ions in solution: Ca40(PO4 )24 (OH)8 (s) + Mg 2 +(aq) → Ca39Mg(PO4 )24 (OH)8 (s) + Ca 2 +(aq)

(6)

The energetic cost of the ion exchange [for the direct reaction in the above equation to the most stable configuration of Mg in the Ca(II) position] is ε = 55.7 kJ/mol. The high value found for ε already indicates a strong energetic preference for the Mg ion to be in water instead of in the HA bulk. We can now estimate the concentration of Mg in HA in equilibrium with aqueous solution as a function of the activities a[Mg2+](aq) and a[Ca2+](aq) of the cations in solution, as proposed in35 xeq =

a[Mg 2 +](aq) a[Ca 2 +](aq)

⎛ ε ⎞ ⎟ exp⎜ − ⎝ RT ⎠

(7) −10

The exponential term at room temperature is 2.0 × 10 , which means that at any reasonable ratio between the concentrations of Mg2+ and Ca2+ in solution, only a very small amount of Mg can be expected to be incorporated in the bulk, assuming equilibrium between the (Mg,Ca)HA solid solution and the aqueous solution. 3.4. Magnesium Substitutions in the Hydrated Surfaces of Hydroxyapatite. The results of the previous sections confirm that the incorporation of Mg into the hydroxyapatite bulk is very unfavorable compared to Mg remaining in the aqueous solution or being located in alternative mineral phases. However, as Mg is found in HA, we now extend our simulations to model the incorporation of Mg into two of the major HA surfaces, (0001) and (011̅0), which are important in the geological36 and biological morphology.37 The PO4−Ca−PO4−PO4−Ca−Ca-terminated (0001) surface is a dipolar surface that can be stabilized by a reconstruction involving the formation of ion vacancies, as shown in Figure 3a. Previous work has shown that this mixed calcium and phosphate termination is marginally more stable than the Ca−Ca termination,38 but, in addition, the PO4−Ca− PO4−PO4−Ca−Ca has a variation of types of cation sites and hence will provide more insight into the different modes of substitutions. There are three surface Ca(II) sites and two types of surface Ca(I) sites in the surface, which differ slightly in their distances to the neighboring oxygen atoms. The (011̅0) plane has one type of cation site [Ca(II)] at the surface, which is coordinated to two oxygen anions. In this surface, the disorder of the hydroxy groups is significant after relaxation of the surface. They remain more or less parallel to the surface plane, although half of the hydroxy groups rotate to interact more closely with neighboring surface calcium ions at Ca(II)OH bond length of 2.4 Å. In the (0001) surface, we observe that the phosphate groups remain intact, although they rotate in order to accommodate the calcium ions. The column of hydroxy groups perpendicular to the surface has become distorted with the topmost group relaxing out of the surface. In the (011̅0) surface, only calcium ions are at the surface; subsurface hydroxy groups are oriented

Figure 3. Segregation energy of magnesium at the surfaces of hydroxyapatite. Ca(I) ions are represented by large green balls in the surface model and by ● symbols in the segregation energy plot, while Ca(II) are represented by small green balls and ■ symbols. The simulation included water molecules, which are not shown here.

parallel to the surface plane, and the disorder is more significant than in the (0001) surface. Futhermore, the surface area of the 011̅0 surface is smaller and will therefore allow less relaxation. We have introduced a single magnesium ion consecutively into all the cation layers of each surface system, also including explicitly in the simulation the water molecules in contact with the solid surfaces. Ca ions were substituted by Mg in the different cation sites of only one side of the slab, but since we are replacing a +2 cation with another +2 cation, this asymmetric slab substitution does not lead to any issues related to surface dipoles. For both surfaces, the magnesium cation showed an overall preference to be incorporated in the Ca(II) sites only in the positions closest to the surface, and this is wellillustrated in Figure 3, which shows the segregation energies versus distance from the surface. The segregation energy ΔEseg is defined as the difference between the energy of Mg incorporated in the surface and the energy of Mg incorporated in the same site in the bulk. Close to the surface the segregation energies are very negative, which means that the magnesium will preferentially locate at the surface rather than in the bulk. Furthermore, these highly negative segregation energies mean that the exchange of Ca for Mg between the surface and aqueous solution (as given by an equation analogous to reaction 6) is exothermic (ε < 0), despite the higher stability of Mg2+(aq) ions with respect to Ca2+(aq) ions in solution. The ion exchange energy for the surface is given by ε(surf) = ε(bulk) + ΔEseg, which leads to ε(0001) = −16.2 kJ/mol and ε(011̅0) = −197.0 kJ/mol. The incorporation of Mg in the HA surfaces is 5854

dx.doi.org/10.1021/la400422d | Langmuir 2013, 29, 5851−5856

Langmuir

Article

clear experimental evidence of preferential incorporation of Mg at the HA surface compared to the bulk.

therefore energetically favorable with respect to ion exchange with aqueous solution, especially in the (011̅0) surface. The distribution of interatomic distances around the cation sites is different for the different surface sites, leading to Mg being more stable when the local changes in MO distances around the cations seem most significant. In the most stable configurations, the coordination distances around Mg in the Ca(II) surface site show a reduction of ∼0.3 Å for Mg−OPO3 and Mg−OH distances compared to Ca−OPO3 and Ca−OH. Hence, the presence of Mg induces changes in the local environment of the different ions, in agreement with experimental observations.9,39 These changes in the local environment are more easily accommodated in the surface than in the bulk material. Moreover, the preference of Mg ions for the surface/water interface can be explained further by their stronger affinity for water compared with Ca. The main limitation of our simulations of the mineral−water interface is the fact that any dissociation of the water molecules, leading to hydroxylation of Ca2+/Mg2+ on the surface and protonation of phosphate groups, is not described in our model. However, ab initio calculations of the interaction of water with MgO and CaO surfaces have shown that the identity of the cation has no significant effect on the reactivity of the surfaces,40−42 but that the dissociation of water at a particular surface site is far more dependent on the geometry, and especially the coordination, of this surface site rather than the identity of the Mg or Ca cation.40,42 As such, our calculations, which are very much a comparison between Ca and Mg sites, are not affected by the mode of adsorption (e.g., associate or dissociative) of the water molecules.

■

AUTHOR INFORMATION

Corresponding Author

*N.A.-B.: e-mail, [email protected]. N.H.D.L.: e-mail, n.h. [email protected]. Present Address ‡

Now at the Institute of Chemical Research of Catalonia ̈ Catalans 16, 43007 Tarragona, Spain. (ICIQ), Av. Paisos Notes

The authors declare no competing financial interest.

■

REFERENCES

(1) Barralet, J.; Best, S.; Bonfield, W. Carbonate substitution in precipitated hydroxyapatite: An investigation into the effects of reaction temperature and bicarbonate ion concentration. J. Biomed. Mater. Res. 1998, 41, 79. (2) Zapantal, R. Effect of carbonate on lattice parameters of apatite. Nature 1965, 206, 403. (3) Jha, L. J.; Best, S. M.; Knowles, J. C.; Rehman, I.; Santos, J. D.; Bonfield, W. Preparation and characterization of fluoride-substituted apatites. J. Mater. Sci: Mater. Med. 1997, 8, 185. (4) Gibson, I. R.; Best, S. M.; Bonfield, W. Chemical characterization of silicon-substituted hydroxyapatite. J. Biomed. Mater. Res. 1999, 44, 422. (5) Bertinetti, L.; Drouet, C.; Combes, C.; Rey, C.; Tampieri, A.; Coluccia, S.; Martra, G. Surface characteristics of nanocrystalline apatites: Effect of Mg surface enrichment on morphology, surface hydration species, and cationic environments. Langmuir 2009, 25, 5647. (6) Drouet, C.; Carayon, M. T.; Combes, C.; Rey, C. Surface enrichment of biomimetic apatites with biologically-active ions Mg2+ and Sr2+: A preamble to the activation of bone repair materials. Mater. Sci. Eng., C 2008, 28, 1544. (7) Dorozhkin, S. V. Calcium orthophosphates. Magnesium deficiency: Effect on bone and mineral metabolism in the mouse. J. Mater. Sci. 2007, 42, 1061. (8) Rude, R. K.; Gruber, H. E.; Wei, L. Y.; Frausto, A.; Mills, B. G. Calcif. Tissue Int. 2003, 72, 32. (9) Laurencin, D.; Almora-Barrios, N.; de Leeuw, N. H.; Gervais, C.; Bonhomme, C.; Mauri, F.; Chrzanowski, W.; Knowles, J. C.; Newport, R. J.; Wong, A.; Gan, Z. H.; Smith, M. E. Magnesium incorporation into hydroxyapatite. Biomaterials 2011, 32, 1826. (10) Bigi, A.; Falini, G.; Foresti, E.; Gazzano, M.; Ripamonti, A.; Roveri, N. Rietveld structure refinements of calcium hydroxylapatite containing magnesium. Acta Crystallogr., Sect. B: Struct. Sci. 1996, 52, 87. (11) Tampieri, A.; Celotti, G.; Landi, E.; M., S. Magnesium doped hydroxyapatite: Synthesis and characterization. Key Eng. Mater. 2004, 264−268, 2051. (12) Sprio, S.; Pezzotti, G.; Celotti, G.; Landi, E.; Tampieri, A. Raman and cathodoluminescence spectroscopies of magnesiumsubstituted hydroxyapatite powders. J. Mater. Res. 2005, 20, 1009. (13) Yasukawa, A.; Ouchi, S.; Kandori, K.; Ishikawa, T. Preparation and characterization of magnesium-calcium hydroxyapatites. J. Mater. Chem. 1996, 6, 1401. (14) Saenger, A. T.; Kuhs, W. F. Structural disorder in hydroxyapatite. Z. Kristallogr. 1992, 199, 123. (15) Stork, L.; Muller, P.; Dronskowski, R.; Ortlepp, J. R. Chemical analyses and X-ray diffraction investigations of human hydroxyapatite minerals from aortic valve stenoses. Z. Kristallogr. 2005, 220, 201. (16) de Leeuw, N. H. Local ordering of hydroxy groups in hydroxyapatite. Chem Commun. (Cambridge, U.K.) 2001, 1646. (17) Corno, M.; Rimola, A.; Bolis, V.; Ugliengo, P. Hydroxyapatite as a key biomaterial: Quantum-mechanical simulation of its surfaces in interaction with biomolecules. Phys. Chem. Chem. Phys. 2010, 12, 6309.

4. CONCLUSIONS We have presented a theoretical investigation of the thermodynamics of Mg incorporation in the bulk and surfaces of hydroxyapatite. While our results show that Mg substitution is energetically more favorable in the Ca(II) than the Ca(I) sites of hydroxyapatite, they also show that Mg incorporation in the HA bulk is very unstable with respect to the formation of separate Mg-rich phases. Mg location in bulk positions is also energetically unfavorable with respect to exchange with ions in aqueous solutions, provided that the Mg and Ca concentrations in solutions are of a similar order of magnitude. Since there is experimental evidence that some Mg does incorporate in the hydroxyapatite bulk, we may conclude that such incorporation is highly metastable and dependent on the growth kinetics, and therefore the cation distribution cannot be expected to be in equilibrium. This rationale can explain why conflicting experimental results have been obtained about the location of Mg in the crystal. Of course, it is possible that the interaction with water molecules present as defects in the hydroxyapatite lattice, which have been identified experimentally by 1H MAS NMR9,43 and the formation of which is energetically only marginally endothermic,44 would cause Mg to be anchored more firmly in the bulk material. However, we have shown that magnesium substitutes calcium much more easily at the hydroxyapatite surfaces, in particular at the (011̅0) termination, where both the energy of cation exchange from solution and the segregation energies are very negative. Exposure of the surface cation sites to the interaction with water favors Mg occupation as Mg is solvated by water more strongly than Ca. Recent IR experiments have reported the presence of Mg ions at the hydroxyapatite surfaces,39 although thus far there is no 5855

dx.doi.org/10.1021/la400422d | Langmuir 2013, 29, 5851−5856

Langmuir

Article

coordinated surface sites of calcium oxide. Phys. Rev. B 2001, 63, 195417. (43) Laurencin, D.; Gervais, C.; Stork, H.; Krämer, S.; Massiot, D.; Fayon, F. 25Mg solid-state NMR of magnesium phosphates: High magnetic field experiments and density functional theory calculations. J. Phys. Chem. C 2012, 116, 19984. (44) de Leeuw, N. H.; Bowe, J. R.; Rabone, J. A. L. A computational investigation of stoichiometric and calcium-deficient oxy- and hydroxyapatite. Faraday Discuss. 2007, 134, 195.

(18) Gale, J. D. GULP: A computer program for the symmetryadapted simulation of solids. J. Chem. Soc., Faraday Trans. 1997, 93, 629. (19) Gale, J. D.; Rohl, A. L. The General Utility Lattice Program (GULP). Mol. Simul. 2003, 29, 291. (20) Gale, J. D. GULP: Capabilities and prospects. Z. Kristallogr. 2005, 220, 552. (21) Born, M.; Huang, K. Dynamical Theory of Crystal Lattices; Clarendon Press: Oxford, 1954. (22) de Leeuw, N. H. Resisting the onset of hydroxyapatite dissolution through the incorporation of fluoride. J. Phys. Chem. B 2004, 108, 1809. (23) Dick, B. G.; Overhauser, A. W. Theory of the dielectric constants of alkali halide crystals. Phys. Rev. 1958, 112, 90. (24) Grau-Crespo, R.; Hamad, S.; Catlow, C. R. A.; de Leeuw, N. H. Symmetry-adapted configurational modelling of fractional site occupancy in solids. J. Phys.-Condens. Matter 2007, 19. (25) Grau-Crespo, R.; de Leeuw, N. H.; Catlow, C. R. A. Cation distribution and magnetic ordering in FeSbO4. J. Mater. Chem. 2003, 13, 2848. (26) Smith, W.; Forester, T. R. DL_POLY_2.0: A general-purpose parallel molecular dynamics simulation package. J. Mol. Graphics 1996, 14, 136. (27) de Leeuw, N. H.; Parker, S. C. Molecular-dynamics simulation of MgO surfaces in liquid water using a shell-model potential for water. Phys. Rev. B: Condens. Matter Mater. Phys. 1998, 58, 13901. (28) de Leeuw, N. H. A computer modelling study of the uptake and segregation of fluoride ions at the hydrated hydroxyapatite (0001) surface: Introducing Ca10(PO4)6(OH)2 potential model. Phys. Chem. Chem. Phys. 2004, 6, 1860. (29) Mkhonto, D.; de Leeuw, N. H. A computer modelling study of the effect of water on the surface structure and morphology of fluorapatite: Introducing a Ca10(PO4)6F2 potential model. J. Mater. Chem. 2002, 12, 2633. (30) Nose, S. Unified formulation of the constant temperature molecular-dynamics methods. J. Chem. Phys. 1984, 81, 511. (31) Hoover, W. G. Canonical dynamics - equilibrium phase - space distributions. Phys. Rev. A 1985, 31, 1695. (32) Verlet, L. Computer experiments on classical fluids I Termodynamical properties of Lennard-Jones molecules. Phys. Rev. 1967, 159, 98. (33) Shannon, R. Revised effective ionic radii and systematic studies of interatomic distances in halides and chalcogenides. Acta Crystallogr., Sect. A 1976, 32, 751. (34) Cotton, F. A.; Wilkinson, G. Advanced Inorganic Chemistry: A Comprehensive Text, 3rd ed.; Interscience: New York, 1972. (35) Ruiz-Hernandez, S. E.; Grau-Crespo, R.; Almora-Barrios, N.; Ruiz-Salvador, A. R.; Fernandez, N.; De Leeuw, N. H. Mg/Ca partitioning between aqueous solution and aragonite mineral: A molecular dynamics study. Chem.−Eur. J. 2012, 18, 9828. (36) Deer, W. A.; Howie, R. A.; Zussman, J. An Introduction to the Rock-Forming Minerals; 2nd ed.; Pearson: Essex, England, 1992. (37) Rohanizadeh, R.; Trecant-Viana, M.; Daculsi, G. Ultrastructural study of apatite precipitation in implanted calcium phosphate ceramic: Influence of the implantation site. Calcif. Tissue Int. 1999, 64, 430. (38) Filgueiras, M. R. T.; Mkhonto, D.; de Leeuw, N. H. Computer simulations of the adsorption of citric acid at hydroxyapatite surfaces. J. Cryst. Growth 2006, 294, 60. (39) Diallo-Garcia, S.; Laurencin, D.; Krafft, J. M.; Casale, S.; Smith, M. E.; Lauron-Pernot, H.; Costentin, G. Influence of magnesium substitution on the basic properties of hydroxyapatites. J. Phys. Chem. C 2011, 115, 24317. (40) Refson, K.; Wogelius, R. A.; Fraser, D. G.; Payne, M. C.; Lee, M. H.; Milman, V. Water chemisorption and reconstruction of the MgO surface. Phys. Rev. B 1995, 52, 10823. (41) Langel, W.; Parrinello, M. Hydrolysis at stepped MgO surfaces. Phys. Rev. Lett. 1994, 73, 504. (42) de Leeuw, N. H.; Purton, J. A. Density-functional theory calculations of the interaction of protons and water with low5856

dx.doi.org/10.1021/la400422d | Langmuir 2013, 29, 5851−5856