**2.1 Defects in hydroxyapatite bulk and surfaces**

Hydroxyapatite (HA) is a mineral which occurs in nature in two polymorphs, a monoclinic form, thermodynamically stable at low temperatures, and an hexagonal form, which can be easily stabilized by substitution of the OH- ions (Suda et al., 1995). These ions are aligned along the c axis (the [001] direction), as highlighted in Fig. 1. The single crystal structure of the hexagonal form of HA is characterized by the *P63/m* space group. The mirror plane, perpendicular to the [001] direction, is compatible with the column of OH ions because of an intrinsic static disorder of these ions, which can point, with no preference, in one of the two opposite directions ([001] or [00-1]). The result is a fractional occupation of the sites in the solved crystallographic structure (50% probability for each direction). As *ab initio* simulation cannot take into account the structural disorder, we reduced the symmetry to *P63*, removing the mirror plane and fixing the directions of the OH ions. In the most stable configuration found, both the OH- ions point in the same direction, as reported in Fig. 1. The oxygen atom of the OH- ion is close to three Ca ions, which form an equilateral triangle in the *ab* plane. Moreover, there are six phosphate ions inside the crystallographic cell, all symmetry equivalent.

The bulk structure of crystalline HA, fully characterized in the literature (Corno et al., 2006), has been considered as a starting point to model the surfaces which are experimentally found to be the most important: (001) in terms of reactivity, and (010) in terms of exposure in the crystal habit (Wierzbicki & Cheung, 2000). Those surfaces have already been fully characterized at an *ab initio* level, and all the structural, geometrical and electronic properties

*In Silico* Study of Hydroxyapatite and Bioglass®: How Computational Science Sheds Light on Biomaterials

Ca/P ratio of 1.62.

using the classical formula:

hydrogen in grey and phosphorous in yellow.

surface as its Ca/P ratio increases to 1.71. The last one is called *P-rich* (010) surface, with a

Fig. 2. The three possible terminations along the [010] direction of HA. Colour coding:

Clearly, due to the non-stoichiometric nature of these new surfaces, Esurf cannot be defined

Esurf = (Eslab - n·Ebulk)/(2·A) (1)

Fig. 3. Upper views of the three (010) surfaces: on the left side, the stoichiometric surface, in the middle, the Ca-rich surface and, on the right side, the P-rich surface. Each structure is superimposed on an isodensity surface colorcoded with the electrostatic potential (blue is most positive, red is most negative). The maps have been generated with the VMD program (Humphrey, Dalke & Schulten, et al., 1996). Calcium represented in cyan, oxygen in red,

calcium cyan, oxygen red, hydrogen grey and phosphorous yellow.

279

have been predicted and compared with the experimental values (Corno et al., 2007). More recently, these modeled surfaces were also employed to study the adsorption of different typologies of molecules that span from water (Corno et al., 2009) to organic acids (Canepa et al., 2011a), amino-acids (Rimola et al., 2008) and small peptides (Corno et al., 2010; Rimola et al., 2009), providing results comparable with the experimental heats of adsorption, when available.

Fig. 1. Two views of the unit cell of the HA bulk structure: calcium in light green, oxygen red, hydrogen light grey and phosphate ions as orange tetrahedra. In the blue cylinders, the column of OH- ions can be observed. Unit cell borders are drawn in light orange.

Nonetheless, the study of the possible terminations along the [010] direction is not yet complete, because two new kinds of (010) surfaces have been discovered experimentally. These surfaces are non-stoichiometric, because their composition is different from the bulk.

#### **2.1.1 The (010) non-stoichiometric surfaces**

In 2002, a HRTEM study highlighted three possible terminations for the (010) surface. As a matter of fact, HA shows an alternation of two layers, Ca3(PO4)2 (A-type) and Ca4(PO4)2(OH)2 (B-type), which can be interrupted in three different ways (Sato et al., 2002). As the sequence is …-A-A-B-A-A-B-A-A-B-…, the periodicity can be truncated by exposing as last layers …-A-B-A or …-A-A-B or …-B-A-A, as shown in Fig. 2. As already done in previous work with other surfaces, we investigated these three possible structures with a slab approach. From the bulk, we extracted a piece of matter along the desired direction, in this case the [010], generating two faces which are exposed to vacuum. These faces ought to be the same to remove any possible dipole moment across the slab, due to geometrical dissimilarities. This necessity may bring the loss of stoichiometricity, in terms of the possibility of obtaining the bulk by replicating the slab along the non-periodic direction: neither B-A-A-B-A-A-B nor A-A-B-A-A-B-A-A can be replicated to regenerate the bulk. Only the stoichiometric surface, characterized by a Ca/P ratio of 1.67, maintains the bulk sequence A-B-A-A-B-A. The slab structure terminating in –A-A-B is called *Ca-rich* (010)

have been predicted and compared with the experimental values (Corno et al., 2007). More recently, these modeled surfaces were also employed to study the adsorption of different typologies of molecules that span from water (Corno et al., 2009) to organic acids (Canepa et al., 2011a), amino-acids (Rimola et al., 2008) and small peptides (Corno et al., 2010; Rimola et al., 2009), providing results comparable with the experimental heats of adsorption, when

Fig. 1. Two views of the unit cell of the HA bulk structure: calcium in light green, oxygen red, hydrogen light grey and phosphate ions as orange tetrahedra. In the blue cylinders, the

Nonetheless, the study of the possible terminations along the [010] direction is not yet complete, because two new kinds of (010) surfaces have been discovered experimentally. These surfaces are non-stoichiometric, because their composition is different from the bulk.

In 2002, a HRTEM study highlighted three possible terminations for the (010) surface. As a matter of fact, HA shows an alternation of two layers, Ca3(PO4)2 (A-type) and Ca4(PO4)2(OH)2 (B-type), which can be interrupted in three different ways (Sato et al., 2002). As the sequence is …-A-A-B-A-A-B-A-A-B-…, the periodicity can be truncated by exposing as last layers …-A-B-A or …-A-A-B or …-B-A-A, as shown in Fig. 2. As already done in previous work with other surfaces, we investigated these three possible structures with a slab approach. From the bulk, we extracted a piece of matter along the desired direction, in this case the [010], generating two faces which are exposed to vacuum. These faces ought to be the same to remove any possible dipole moment across the slab, due to geometrical dissimilarities. This necessity may bring the loss of stoichiometricity, in terms of the possibility of obtaining the bulk by replicating the slab along the non-periodic direction: neither B-A-A-B-A-A-B nor A-A-B-A-A-B-A-A can be replicated to regenerate the bulk. Only the stoichiometric surface, characterized by a Ca/P ratio of 1.67, maintains the bulk sequence A-B-A-A-B-A. The slab structure terminating in –A-A-B is called *Ca-rich* (010)

column of OH- ions can be observed. Unit cell borders are drawn in light orange.

**2.1.1 The (010) non-stoichiometric surfaces** 

available.

surface as its Ca/P ratio increases to 1.71. The last one is called *P-rich* (010) surface, with a Ca/P ratio of 1.62.

Fig. 2. The three possible terminations along the [010] direction of HA. Colour coding: calcium cyan, oxygen red, hydrogen grey and phosphorous yellow.

Clearly, due to the non-stoichiometric nature of these new surfaces, Esurf cannot be defined using the classical formula:

$$\mathbf{E}\_{\text{surf}} = (\mathbf{E}\_{\text{slab}} \cdot \mathbf{n} \cdot \mathbf{E}\_{\text{bulk}}) / \text{(2:A)}\tag{1}$$

Fig. 3. Upper views of the three (010) surfaces: on the left side, the stoichiometric surface, in the middle, the Ca-rich surface and, on the right side, the P-rich surface. Each structure is superimposed on an isodensity surface colorcoded with the electrostatic potential (blue is most positive, red is most negative). The maps have been generated with the VMD program (Humphrey, Dalke & Schulten, et al., 1996). Calcium represented in cyan, oxygen in red, hydrogen in grey and phosphorous in yellow.

*In Silico* Study of Hydroxyapatite and Bioglass®: How Computational Science Sheds Light on Biomaterials

adducts are reported for water and CO.

for a loading of water comparable with our models.

symmetric mode and *vice versa* (Corno et al., 2009).

2007).

almost equivalence for the two Ca ions for the two surfaces.

molecules, those commonly used are carbon monoxide and water, because of their selectivity towards cations, the ease of interpretation of their spectra, and the high sensitivity of IR instrumentation in the region where their vibrational frequencies fall (medium IR). The (010) stoichiometric surface has already been fully characterized by Corno *et al.* in relation to the adsorption of these molecules (Corno et al., 2009), whereas the nonstoichiometric surfaces are the subject of this work. Calculations provide binding energies and vibrational frequencies of the probe molecules which can be used as a future reference to be compared with experimental measurements. In Fig. 4, the optimized structures of the

The adsorptions of water upon the most exposed cations of the non-stoichiometric surfaces are characterized by BE values of 131 and 125 kJ/mol (BSSE ≈ 35 %), for the Ca-rich and Prich surfaces, respectively. These values also take into account the formation of two hydrogen bonds between the water molecule and exposed anions, either phosphate only, or also hydroxyl anions in the case of the (010) Ca-rich surface. If the thermal and vibrational contributions are taken into account, these values decrease to 117 and 110 kJ/mol, respectively, allowing a consistent comparison with the experimental differential heat of adsorption of water of 110 kJ/mol (Corno et al., 2009). The latter value has been obtained from experimental micro-calorimetric studies of water adsorption on microcrystalline HA,

The adsorption of CO upon the most exposed Ca ions gives results highly representative of the strength of the cationic site as no hydrogen bond can be formed: the BE values upon the Ca-rich and P-rich surfaces are, respectively, 38 and 40 kJ/mol (BSSE ≈ 20 %), showing an

When the vibrational features are considered, a crucial point is the comparison between the frequencies of the free and the adsorbate molecule. While the stretching mode of CO is easily identified in both cases, the modes of the adsorbed water molecule are no longer easily referable to those of the free molecule. The free water molecule is characterized by two stretching modes (the symmetric and, at higher frequencies, the anti-symmetric) and one bending mode. When the molecule interacts with the surface, these two kinds of stretching modes are no longer classifiable on a symmetry ground. The hydrogen bonds which are formed with the most exposed anions of the surface cause the loss of the C2v symmetry of the molecule: each OH bond oscillates independently. As in previous work, we decided to compare the lower OH stretching frequency of the adsorbed molecule to the

With these remarks, the calculated symmetric stretching shifts are -761 and -329 cm-1 for the (010) Ca-rich and (010) P-rich surfaces, respectively, while the anti-symmetric shifts are -310 and -328 cm-1. The shifts of the bending frequencies are 113 and 88 cm-1. When a hydrogen bond pattern occurs, the stretching mode shifts are negative due to the electronic density transfer and the resulting decrease of the bond strength. Instead, the bending frequencies increase because of the restraint caused by the hydrogen bond itself. The larger shift of the symmetric stretching of the H2O adsorbed on the (010) Ca-rich surface, in combination with a larger bending shift, indicates the formation of a very strong hydrogen bond with an exposed anion of the (010) Ca-rich surface. Experimentally, the shift of the stretching is -400 cm-1 while the bending shift is 40 cm-1. The calculated values are, then, in agreement with the experimental ones as long as the signs and the trends are considered (Bertinetti et al.,

281

where A is the surface area (doubled, because two faces are generated), Ebulk and Eslab are the calculated energies of the related models and n is an integer number, required to match the chemical potentials of the two systems. Astala and Stott adopted a clever and rather involved scheme to evaluate the phase existence conditions for the two non-stoichiometric surfaces, showing that the region of their stability is outside the stability window defined by bulk HA, Ca(OH)2 and -Ca3(PO4)2 (Astala & Stott, 2008).

The typology of the exposed layers differentiates the (010) stoichiometric surface from the non-stoichiometric ones. Indeed, the (010) Ca-rich surface exposes the OH ions belonging to the B layer, while both the stoichiometric and the P-rich surfaces contain the OH ions in an inner layer (see the structures displayed in Fig. 3 for details).


\* The slab thickness is the perpendicular distance between the most exposed Ca ions on the upper and lower faces of the slab.

Table 1. Most important geometrical parameters of the HA (010) surfaces. The values in <…> correspond to the arithmetical averages of the considered feature.

In Fig. 3, the most exposed layers of the three terminations are reported, superimposed on the isodensity surface colorcoded with the electrostatic potential. The isodensity value is fixed to 10-6 electrons and the electrostatic potential spans from -0.02 a.u. (red) to +0.02 a.u. (blue). Positive values of the potential are visible in correspondence of the Ca ions while negative potential zones are mostly located upon superficial phosphate ions.

In Table 1, a geometrical analysis of the three surfaces is reported. The cell parameters are slightly different between the three cases, but the rectangular shape is mostly maintained. The slab thickness is not comparable, because of a different number of layers for each surface model. The interatomic distances and angles are, however, very similar.

Two important intensive parameters can classify the stability of a slab structure, the electronic band gap and the total dipole moment across the slab: each surface has a dipole moment close to zero and a band gap typical of electrical insulators. The three surfaces are, then, stable and can be adopted as a substrate to study the adsorption of molecules.

#### **2.1.2 Adsorption of H2O and CO upon the most exposed Ca ions of the (010) surfaces**

The chemical reactivity of the most exposed cations of a surface is experimentally studied with the IR technique by monitoring the perturbation of the vibrational frequency of the probe molecule, which is compared to the value for the free molecule. Of the possible probe

where A is the surface area (doubled, because two faces are generated), Ebulk and Eslab are the calculated energies of the related models and n is an integer number, required to match the chemical potentials of the two systems. Astala and Stott adopted a clever and rather involved scheme to evaluate the phase existence conditions for the two non-stoichiometric surfaces, showing that the region of their stability is outside the stability window defined by

The typology of the exposed layers differentiates the (010) stoichiometric surface from the

the B layer, while both the stoichiometric and the P-rich surfaces contain the OH ions in an

Table 1. Most important geometrical parameters of the HA (010) surfaces. The values in <…>

In Fig. 3, the most exposed layers of the three terminations are reported, superimposed on the isodensity surface colorcoded with the electrostatic potential. The isodensity value is fixed to 10-6 electrons and the electrostatic potential spans from -0.02 a.u. (red) to +0.02 a.u. (blue). Positive values of the potential are visible in correspondence of the Ca ions while

In Table 1, a geometrical analysis of the three surfaces is reported. The cell parameters are slightly different between the three cases, but the rectangular shape is mostly maintained. The slab thickness is not comparable, because of a different number of layers for each

Two important intensive parameters can classify the stability of a slab structure, the electronic band gap and the total dipole moment across the slab: each surface has a dipole moment close to zero and a band gap typical of electrical insulators. The three surfaces are,

**2.1.2 Adsorption of H2O and CO upon the most exposed Ca ions of the (010) surfaces**  The chemical reactivity of the most exposed cations of a surface is experimentally studied with the IR technique by monitoring the perturbation of the vibrational frequency of the probe molecule, which is compared to the value for the free molecule. Of the possible probe

a (Å) 6.98 6.93 7.00 b (Å) 9.28 9.27 9.33 γ (°) 89.87 90.01 90.02 Area (Å2) 64.83 64.24 65.24 \* Thickness (Å) 13.28 20.05 20.23 <Ca-O> (Å) 2.41 2.39 2.40 <P-O> (Å) 1.55 1.55 1.55 <O-H> (Å) 0.97 0.97 0.97 <OPO> (°) 109.5 109.4 109.4 Band Gap (eV) 6.21 6.80 6.65 Dipole (Debye) 3.5·10-3 4.0·10-3 1.4·10-2 \* The slab thickness is the perpendicular distance between the most exposed Ca ions on

Stoichiometric Ca-rich P-rich

ions belonging to

bulk HA, Ca(OH)2 and -Ca3(PO4)2 (Astala & Stott, 2008).

inner layer (see the structures displayed in Fig. 3 for details).

correspond to the arithmetical averages of the considered feature.

negative potential zones are mostly located upon superficial phosphate ions.

surface model. The interatomic distances and angles are, however, very similar.

then, stable and can be adopted as a substrate to study the adsorption of molecules.

the upper and lower faces of the slab.

non-stoichiometric ones. Indeed, the (010) Ca-rich surface exposes the OH-

molecules, those commonly used are carbon monoxide and water, because of their selectivity towards cations, the ease of interpretation of their spectra, and the high sensitivity of IR instrumentation in the region where their vibrational frequencies fall (medium IR). The (010) stoichiometric surface has already been fully characterized by Corno *et al.* in relation to the adsorption of these molecules (Corno et al., 2009), whereas the nonstoichiometric surfaces are the subject of this work. Calculations provide binding energies and vibrational frequencies of the probe molecules which can be used as a future reference to be compared with experimental measurements. In Fig. 4, the optimized structures of the adducts are reported for water and CO.

The adsorptions of water upon the most exposed cations of the non-stoichiometric surfaces are characterized by BE values of 131 and 125 kJ/mol (BSSE ≈ 35 %), for the Ca-rich and Prich surfaces, respectively. These values also take into account the formation of two hydrogen bonds between the water molecule and exposed anions, either phosphate only, or also hydroxyl anions in the case of the (010) Ca-rich surface. If the thermal and vibrational contributions are taken into account, these values decrease to 117 and 110 kJ/mol, respectively, allowing a consistent comparison with the experimental differential heat of adsorption of water of 110 kJ/mol (Corno et al., 2009). The latter value has been obtained from experimental micro-calorimetric studies of water adsorption on microcrystalline HA, for a loading of water comparable with our models.

The adsorption of CO upon the most exposed Ca ions gives results highly representative of the strength of the cationic site as no hydrogen bond can be formed: the BE values upon the Ca-rich and P-rich surfaces are, respectively, 38 and 40 kJ/mol (BSSE ≈ 20 %), showing an almost equivalence for the two Ca ions for the two surfaces.

When the vibrational features are considered, a crucial point is the comparison between the frequencies of the free and the adsorbate molecule. While the stretching mode of CO is easily identified in both cases, the modes of the adsorbed water molecule are no longer easily referable to those of the free molecule. The free water molecule is characterized by two stretching modes (the symmetric and, at higher frequencies, the anti-symmetric) and one bending mode. When the molecule interacts with the surface, these two kinds of stretching modes are no longer classifiable on a symmetry ground. The hydrogen bonds which are formed with the most exposed anions of the surface cause the loss of the C2v symmetry of the molecule: each OH bond oscillates independently. As in previous work, we decided to compare the lower OH stretching frequency of the adsorbed molecule to the symmetric mode and *vice versa* (Corno et al., 2009).

With these remarks, the calculated symmetric stretching shifts are -761 and -329 cm-1 for the (010) Ca-rich and (010) P-rich surfaces, respectively, while the anti-symmetric shifts are -310 and -328 cm-1. The shifts of the bending frequencies are 113 and 88 cm-1. When a hydrogen bond pattern occurs, the stretching mode shifts are negative due to the electronic density transfer and the resulting decrease of the bond strength. Instead, the bending frequencies increase because of the restraint caused by the hydrogen bond itself. The larger shift of the symmetric stretching of the H2O adsorbed on the (010) Ca-rich surface, in combination with a larger bending shift, indicates the formation of a very strong hydrogen bond with an exposed anion of the (010) Ca-rich surface. Experimentally, the shift of the stretching is -400 cm-1 while the bending shift is 40 cm-1. The calculated values are, then, in agreement with the experimental ones as long as the signs and the trends are considered (Bertinetti et al., 2007).

*In Silico* Study of Hydroxyapatite and Bioglass®: How Computational Science Sheds Light on Biomaterials

The results of these studies assert that the CO3

substitution.

reaction (2):

compensation is required.

ways (Astala & Stott, 2005):

on the hydroxyapatite structure, in order to comprehend where and how this ion is located.

*type A* defect, or a PO43-, a *type B* defect. This hypothesis has been confirmed by IR spectra in which the vibrational mode frequencies of the carbonate are different for each type of

As the net charge of the carbonate is different from those of the substitutional anions, charge

In the case of the type A defect, the charge excess can be compensated by creation of an ionic vacancy removing another OH ion, also because of the larger steric encumbrance of the

In the case of the type B defect, the electro-neutrality can be maintained in five different

When considering the type A defect, it is first necessary to notice that, in the hydroxyapatite unit cell, only two OH ions are present (Fig. 1). If they are both removed at once but only one carbonate ion is included, the resulting structure is no longer a hydroxyapatite, as it had lost all the hydroxyl ions, and has to be considered a carbonated apatite. In Fig. 5, this stable structure is reported. The similarity between carbonate apatite and calcite structures is clear: the averaged distance <Ca-OC> is 2.35 Å in the carbonate apatite and 2.34 Å in the calcite. From the *ab initio* calculation, it is also possible to predict the enthalpy of formation of the carbonate apatite from the calcium phosphate and the calcite structures, considering the

 3Ca3(PO4)2 + CaCO3 → Ca10(PO4)6(CO3) (2) The B3LYP enthalpy of formation is -35 kJ/mol, directly comparable with the value of -32 kJ/mol obtained with the VASP code and the PW91 functional (Rabone & de Leeuw, 2007).

Fig. 5. Carbonated apatite. On the left side, a view of the carbonate ion is reported, highlighting its distances to the closest Ca ions. On the right side, the cell, where both OH are removed and one carbonate is substituted, is reported after full relaxation of the atoms.

carbonate ion with respect to the hydroxyl ion (Peroos et al., 2006).

3. Substitution of one Ca with one hydrogen (B3 complex);

4. Substitution of one Ca with one alkaline ion; 5. OH ion incorporation close to the carbonate ion.

1. Removal of one OH ion and creation of one Ca vacancy (B1 complex); 2. Removal of one Ca ion for every two carbonate ions (B2 complex);

283

2- can substitute an OH- ion of the structure, a

Fig. 4. Best views of the interactions between the most exposed ions of the (010) Ca-rich and (010) P-rich surfaces and the two probe molecules, water and carbon monoxide. The interactions with Ca ions and the hydrogen bonds are represented with dotted lines.

The stretching frequency shifts of the CO molecule are 40 cm-1 for the (010) Ca-rich surface and 41 cm-1 for the (010) P-rich surface. These BE and frequency shifts values indicate that the Ca ion can reasonably be considered as the major cause of the activity of HA towards polar molecules free from H-bond interactions.

A comparison with the experimental IR spectrum of CO adsorbed on HA reveals two main components, characterized by an upward shift of the stretching frequency of the adsorbed molecule of about +27 and +41 cm-1, respectively. The highest shift is attributed to a tiny fraction of the most exposed calcium ions, which are those modeled in our studies, so proving a good agreement. On the other hand, the majority of calcium sites contribute to the lowest shift (Bertinetti et al., 2007; Sakhno et al., 2010).

#### **2.1.3 Hydroxyapatite and carbonate ion defects**

HA is the main component of the inorganic phase of bones and teeth, but, as many studies have demonstrated, the mineral present in those tissues is neither crystalline nor without defects (Fleet & Liu, 2003; Fleet & Liu, 2004; Fleet & Liu, 2007; Rabone & de Leeuw, 2005; Astala & Stott, 2005; Astala et al., 2006; Rabone & de Leeuw, 2007; de Leeuw et al., 2007 ). Indeed, it incorporates many other elements, ions or compounds, such as Mg or other alkaline earth metals instead of Ca, and, above all, the carbonate anion. Many studies, both theoretical and experimental, have already been conducted on the role of the carbonate ion

Fig. 4. Best views of the interactions between the most exposed ions of the (010) Ca-rich and

The stretching frequency shifts of the CO molecule are 40 cm-1 for the (010) Ca-rich surface and 41 cm-1 for the (010) P-rich surface. These BE and frequency shifts values indicate that the Ca ion can reasonably be considered as the major cause of the activity of HA towards

A comparison with the experimental IR spectrum of CO adsorbed on HA reveals two main components, characterized by an upward shift of the stretching frequency of the adsorbed molecule of about +27 and +41 cm-1, respectively. The highest shift is attributed to a tiny fraction of the most exposed calcium ions, which are those modeled in our studies, so proving a good agreement. On the other hand, the majority of calcium sites contribute to the

HA is the main component of the inorganic phase of bones and teeth, but, as many studies have demonstrated, the mineral present in those tissues is neither crystalline nor without defects (Fleet & Liu, 2003; Fleet & Liu, 2004; Fleet & Liu, 2007; Rabone & de Leeuw, 2005; Astala & Stott, 2005; Astala et al., 2006; Rabone & de Leeuw, 2007; de Leeuw et al., 2007 ). Indeed, it incorporates many other elements, ions or compounds, such as Mg or other alkaline earth metals instead of Ca, and, above all, the carbonate anion. Many studies, both theoretical and experimental, have already been conducted on the role of the carbonate ion

(010) P-rich surfaces and the two probe molecules, water and carbon monoxide. The interactions with Ca ions and the hydrogen bonds are represented with dotted lines.

polar molecules free from H-bond interactions.

lowest shift (Bertinetti et al., 2007; Sakhno et al., 2010).

**2.1.3 Hydroxyapatite and carbonate ion defects** 

on the hydroxyapatite structure, in order to comprehend where and how this ion is located. The results of these studies assert that the CO3 2- can substitute an OH- ion of the structure, a *type A* defect, or a PO43-, a *type B* defect. This hypothesis has been confirmed by IR spectra in which the vibrational mode frequencies of the carbonate are different for each type of substitution.

As the net charge of the carbonate is different from those of the substitutional anions, charge compensation is required.

In the case of the type A defect, the charge excess can be compensated by creation of an ionic vacancy removing another OH ion, also because of the larger steric encumbrance of the carbonate ion with respect to the hydroxyl ion (Peroos et al., 2006).

In the case of the type B defect, the electro-neutrality can be maintained in five different ways (Astala & Stott, 2005):


When considering the type A defect, it is first necessary to notice that, in the hydroxyapatite unit cell, only two OH ions are present (Fig. 1). If they are both removed at once but only one carbonate ion is included, the resulting structure is no longer a hydroxyapatite, as it had lost all the hydroxyl ions, and has to be considered a carbonated apatite. In Fig. 5, this stable structure is reported. The similarity between carbonate apatite and calcite structures is clear: the averaged distance <Ca-OC> is 2.35 Å in the carbonate apatite and 2.34 Å in the calcite. From the *ab initio* calculation, it is also possible to predict the enthalpy of formation of the carbonate apatite from the calcium phosphate and the calcite structures, considering the reaction (2):

$$\text{Ca}\_3\text{(PO}\_4\text{)}\_2 + \text{CaCO}\_3 \rightarrow \text{Ca}\_{10}\text{(PO}\_4\text{)}\_6\text{(CO}\_3\text{)}\tag{2}$$

The B3LYP enthalpy of formation is -35 kJ/mol, directly comparable with the value of -32 kJ/mol obtained with the VASP code and the PW91 functional (Rabone & de Leeuw, 2007).

Fig. 5. Carbonated apatite. On the left side, a view of the carbonate ion is reported, highlighting its distances to the closest Ca ions. On the right side, the cell, where both OH are removed and one carbonate is substituted, is reported after full relaxation of the atoms.

*In Silico* Study of Hydroxyapatite and Bioglass®: How Computational Science Sheds Light on Biomaterials

typical of classical techniques (Tilocca, 2010).

theoretical research work has been carried out on bioglasses using *ab initio* molecular dynamics, though the transferability of empirical potentials remains a challenging goal,

Fig. 6. Top views of the most stable carbonated hydroxyapatite structures. The carbonate can substitute a pair of hydroxyl ions (green circles) or a phosphate ion (blue circles). Top left, type A defect with three out of four pairs of OH ions substituted; top right, type B defect with Na as charge compensator; bottom left, type B with H as charge compensator; bottom

The *ab initio* modelling of an amorphous material as the 45S5 Bioglass® has straightway caused a number of challenges. As the internal coordinates are not available from classical structural analysis, we have developed a multi-scale strategy to obtain a feasible amorphous bioglass model similar in composition to the 45S5. This approach has been largely illustrated in previous papers (Corno & Pedone, 2009; Corno et al., 2008; A. Pedone et al., 2008) so that

Firstly, we adopted classical molecular dynamics simulations to model a glassy bulk structure which could be close to the 45S5 composition keeping the size of the unit cell small enough for *ab initio* calculations, *i.e.* 78 atoms. In order to reproduce the correct ratio between SiO2, P2O5, Na2O and CaO components, the experimental density of 2.72 g/cm3 has been maintained fixed in a cubic box of 10.10 Å per side. After randomly generating atomic

right, type A and B, with H as charge compensator of the B type.

**2.2.1 Multi-scale strategy for the modelling of 45S5 Bioglass®** 

here we only summarize the most important steps.

285

As the main interest in the type A defect is the stability of the OH- -CO3 2- substitution as a function of CO32- content, we modeled three cases using a 2x2 supercell:

$$\text{(4-x)Cau(PO4)\_6(OH)\_2 + xCa\_{10}(PO4)\_6(CO)\_2 \rightarrow Ca\_{40}(PO4)\_{24}(CO)\_2(OH)\_{8\cdot 24} \tag{3}$$

The most stable case occurs (see Fig. 6A) when three out of four pairs of OH ions (x=3) were substituted with carbonate ions (ΔEr = -24 kJ/mol). It is notable that the percentage of carbonate related to this minimum energy structure (4.4%) is close to that found in dentine tissue (5.6%) and in tooth enamel (3.5%) (Dorozhkin, 2009).

Among the five typologies of B type defect, we investigated only the substitution of Ca with Na or H. As all the phosphate ions are equivalent by symmetry, the substitution with carbonate is easy. On the contrary, in the cationic substitution to restore the neutrality, there is a choice among ten Ca ions. These are not equivalent, as the symmetry is broken by the carbonate substitution. We selected only four Ca ions, those nearest to the carbonate, because from the experiment it is known that the two defective substituents are close to each other. Indeed, in the four selected cases, the most stable structures are those in which the distance between the substituents is minimum. The different stabilities were obtained by calculating the energy variation (ΔEr) for the following reactions:

$$\text{Ca}\_{10}\text{(PO}\_4\text{)}\_6\text{(OH)}\_2 + \text{NaHCO}\_3 \rightarrow \text{Ca}\_9\text{Na(PO}\_4\text{)}\_5\text{(CO}\_3\text{)(OH)}\_2 + \text{CaHPO}\_4 \tag{4}$$

$$\text{Ca}\_{10}(\text{PO}\_4)\_6(\text{OH})\_2 + \text{H}\_2\text{CO}\_3 \rightarrow \text{CaAl(PO}\_4)\_5(\text{CO}\_2)(\text{OH})\_2 + \text{CaHPO}\_4 \tag{5}$$

The ΔEr is 108 kJ/mol for reaction (4) and -94 kJ/mol for reaction (5), indicating that the inclusion of H is much more preferable than Na. The reason relies upon the formation of bulk water between the H and the OH ion of the column. This water molecule stabilizes the structure and its removal requires 136 kJ/mol, because of the occurrence of rather strong Hbonds. The two structures are reported in Fig. 6 (B/Na and B/H).

We also calculated some bulk structures in which both typologies (A and B) of substitutions were present at the same time, in order to better mimic the bone features. Among all the simulated structures, the most favorable situation is the one reported in Fig. 6 (A+B/H): a carbonate ion substitutes a pair of OH ions while another carbonate replaces a phosphate with formation of a water molecule. The hydrogen bonds formed by the water molecule are stronger than those of the B defect model, highlighting that the two defects interact and influence each other. The energy of formation of the mixed A+B/H structure is -752 kJ/mol for the reaction 6.

$$11.5\,\mathrm{Ca(PO\_4)\_2} + 1.5\,\mathrm{CaCO\_3} + 3\,\mathrm{Ca(OH)\_2} + 0.5\,\mathrm{H\_2CO\_3} \to \mathrm{Ca9H(PO\_4)\_2(CO\_3)\_8(OH)\_6} \tag{6}$$

#### **2.2 Bioglass: the effect of varying phosphorous content**

As already described in the Introduction of this Chapter, bioactive glasses are extensively studied as prostheses for bone and tooth replacement and regeneration. In particular, the 45S5 composition has continuously been investigated, not only in its compact form, whose applications are limited to low load-bearing, but also as particulates and powders for bone filler use. Hence, the computational investigations of the variation in composition for the different glass components still represent a very interesting and stimulating task.

In the computational area, the most natural method to simulate glassy materials is usually classical molecular dynamics via the melt-and-quench procedure. Recently, much

 (4-x)Ca10(PO4)6(OH)2 + xCa10(PO4)6(CO3) → Ca40(PO4)24(CO3)x(OH)8-2x (3) The most stable case occurs (see Fig. 6A) when three out of four pairs of OH ions (x=3) were substituted with carbonate ions (ΔEr = -24 kJ/mol). It is notable that the percentage of carbonate related to this minimum energy structure (4.4%) is close to that found in dentine

Among the five typologies of B type defect, we investigated only the substitution of Ca with Na or H. As all the phosphate ions are equivalent by symmetry, the substitution with carbonate is easy. On the contrary, in the cationic substitution to restore the neutrality, there is a choice among ten Ca ions. These are not equivalent, as the symmetry is broken by the carbonate substitution. We selected only four Ca ions, those nearest to the carbonate, because from the experiment it is known that the two defective substituents are close to each other. Indeed, in the four selected cases, the most stable structures are those in which the distance between the substituents is minimum. The different stabilities were obtained by

Ca10(PO4)6(OH)2 + NaHCO3 → Ca9Na(PO4)5(CO3)(OH)2 + CaHPO4 (4)

 Ca10(PO4)6(OH)2 + H2CO3 → Ca9H(PO4)5(CO3)(OH)2 + CaHPO4 (5) The ΔEr is 108 kJ/mol for reaction (4) and -94 kJ/mol for reaction (5), indicating that the inclusion of H is much more preferable than Na. The reason relies upon the formation of bulk water between the H and the OH ion of the column. This water molecule stabilizes the structure and its removal requires 136 kJ/mol, because of the occurrence of rather strong H-

We also calculated some bulk structures in which both typologies (A and B) of substitutions were present at the same time, in order to better mimic the bone features. Among all the simulated structures, the most favorable situation is the one reported in Fig. 6 (A+B/H): a carbonate ion substitutes a pair of OH ions while another carbonate replaces a phosphate with formation of a water molecule. The hydrogen bonds formed by the water molecule are stronger than those of the B defect model, highlighting that the two defects interact and influence each other. The energy of formation of the mixed A+B/H structure is -752 kJ/mol

11.5 Ca3(PO4)2 + 1.5 CaCO3 + 3 Ca(OH)2 + 0.5 H2CO3 → Ca39H(PO4)23(CO3)2(OH)6 (6)

As already described in the Introduction of this Chapter, bioactive glasses are extensively studied as prostheses for bone and tooth replacement and regeneration. In particular, the 45S5 composition has continuously been investigated, not only in its compact form, whose applications are limited to low load-bearing, but also as particulates and powders for bone filler use. Hence, the computational investigations of the variation in composition for the

In the computational area, the most natural method to simulate glassy materials is usually classical molecular dynamics via the melt-and-quench procedure. Recently, much

different glass components still represent a very interesting and stimulating task.


As the main interest in the type A defect is the stability of the OH-

tissue (5.6%) and in tooth enamel (3.5%) (Dorozhkin, 2009).

calculating the energy variation (ΔEr) for the following reactions:

bonds. The two structures are reported in Fig. 6 (B/Na and B/H).

**2.2 Bioglass: the effect of varying phosphorous content** 

for the reaction 6.

function of CO32- content, we modeled three cases using a 2x2 supercell:

theoretical research work has been carried out on bioglasses using *ab initio* molecular dynamics, though the transferability of empirical potentials remains a challenging goal, typical of classical techniques (Tilocca, 2010).

Fig. 6. Top views of the most stable carbonated hydroxyapatite structures. The carbonate can substitute a pair of hydroxyl ions (green circles) or a phosphate ion (blue circles). Top left, type A defect with three out of four pairs of OH ions substituted; top right, type B defect with Na as charge compensator; bottom left, type B with H as charge compensator; bottom right, type A and B, with H as charge compensator of the B type.
