Atomistic Simulation of Intrinsic Defects and Trivalent and Tetravalent Ion Doping in Hydroxyapatite

Atomistic simulation techniques have been employed in order to investigate key issues related to intrinsic defects and a variety of dopants from trivalent and tetravalent ions. The most favorable intrinsic defect is determined to be a scheme involving calcium and hydroxyl vacancies. It is found that trivalent ions have an energetic preference for the Ca site, while tetravalent ions can enter P sites. Charge compensation is predicted to occur basically via three schemes. In general, the charge compensation via the formation of calcium vacancies is more favorable. Trivalent dopant ions are more stable than tetravalent dopants.


Introduction
Hydroxyapatite (Ca 10 (PO 4 ) 6 (OH) 2 ) has been intensively investigated due to the possibility of applying it in a range of biomedical applications.This is attributed to characteristics such as its mechanical properties, biocompatibility, osteoconductivity, nontoxic properties, noninflammatory properties, and the similarity of the mineral constituents to human bone and teeth [1][2][3][4][5][6].Due to these properties, it can be applied in the manufacture of artificial bone material, as a coating on surgical implants, and in tissue engineering, drug and gene delivery, and other biomedical areas [7][8][9][10][11].Hydroxyapatite has a hexagonal crystal structure, with space group P63/m, where the experimental dimensions of the unit cell are  =  = 9.419 Å,  = 6.881Å,  =  = 90 ∘ , and  = 120 ∘ [12].There are two crystallographically different sites for Ca 2+ , which have different coordination numbers.Ca(I) is coordinated with nine oxygen ions and Ca(II) is coordinated with six oxygen ions and one from the hydroxyl group.The average Ca-O distances are lesser than 3 Å.
Many works have investigated the substitution of extrinsic dopant in hydroxyapatite (HAP) in order to improve the material properties [13][14][15][16][17].The incorporation of dopant in the HAP structure contributes to increasing a range of potential applications in fields such as water purification, bone pathologies, bioceramics, catalysis, and luminescence [18].For example, it is reported that Sb 3+ -doped HAP has useful applications in fluorescent lamps [19].Cr 3+ was incorporated as a dopant in the material with the aim of developing a biosensor [20].HAP can also be doped by trivalent lanthanide ions such as Nd 3+ [21,22], Yb 3+ [23,24], Er 3+ [25], Gd 3+ , and Eu 3+ [26,27].The incorporation of Eu and Gd in a HAP structure can affect the bone remodeling cycle, and they have potential for the treatment of bone density disorders such as osteoporosis [26,27].Eu ions exhibit favorable spectroscopic properties and can also be used for biological imaging [28,29] and in lasers.On the other hand, gadolinium can be used as a contrast agent to provide a brighter magnetic resonance (MR) imaging signal [30].
A systematic investigation of cation substitution in HAP may provide a better understanding of their properties, which is still lacking.Some theoretical studies based on first principle calculations for HAP materials have been realized in order to elucidate some basic physics [31][32][33].However, many questions regarding the study of defects remain open.Dopant ions with 3+ and 4+ valence states such as rare earth ions and some transition metals will be incorporated into the Ca or P sites.The compensation mechanism for this doping remains inconclusive.For this reason, this paper makes use of computer modeling methods in order to predict the preferred dopant sites and charge compensation mechanisms in the material.
For Eu-doped HAP, for example, an investigation of the structural response is essential due to the intimate relationship between its properties and lattice distortions caused by the Eu substitution and the charge compensation involved in the Eu emission, which are associated with the crystal field [34].It is well known that the emission spectrum of Eu 3+ strongly depends on the symmetry of the site which Eu 3+ occupies.If it occupies a site with inversion symmetry, only the magnetic-dipole transition 5 D 0 -7 F 1 can be observed; but if there is no inversion symmetry at the site of the Eu 3+ ion, the electric-dipole transitions 5 D 0 -7 F 2 can be observed.Also, the transition is sensitive to the coordination environment [35].
In the present work, atomistic simulation based on the lattice energy minimization is used to provide useful information on the energetically favoured intrinsic defects and a variety of dopants from trivalent and tetravalent ions.Specifically, trivalent rare earth ions and some trivalent and tetravalent transition metals were considered due to their influence on various properties [31][32][33].A detailed investigation on the defects behavior of HAP is very important in order to understand their importance on the improvement of the material properties.In these methods, detailed estimation of lattice relaxation and the Coulomb energies of a large number of lattice ions around defect species are provided without high computational demanding, on the contrary of the DFT based methods.

Methodology
The simulation techniques used in this work are based on energy minimization, with interactions represented by interatomic potentials.Interactions between ions in the solid are represented in terms of a long-range Coulomb term plus a short-range term, as described by Buckingham potential, that accounts for electron cloud overlap Pauli repulsion and dispersion (Van der Waals) interactions.Polarizability of the oxygen ions is incorporated by means of the Dick-Overhauser shell model [36].Most of the interatomic potential models that were obtained empirically for the simulation of hydroxyapatite materials [37][38][39][40][41] showed excellent agreement with the experimental structures and properties.However, in order to improve the fit with respect to structural and mechanical properties, the interatomic potential initially derived by Mostafa and Brown [42] was refitted empirically using the GULP code [43].In our refit, the  and  parameters for P-O interaction from Mostafa and Brown [42] were modified.In our work, the tree body potential for P-O interaction used by Mostafa and Brown [42] was not considered.Defects are modeled using the Mott-Littleton approximation [44], in which a spherical region of lattice surrounding the defect is treated explicitly, with all interactions being considered, and more distant parts of the lattice are treated using a continuum approach.

Results
Suitable potential parameters are highly necessary in order to give a good description of the material by atomistic simulation.The refitted potential parameters used in this work are listed in Table 1.The obtained structural parameters are compared with the experimental values and other theoretical calculations in Table 2.It can be seen that the lattice parameters obtained in this work are more accurate than those obtained previously using atomistic simulation with the refit procedure [42] or using other force fields with the partial charge model [45][46][47] or those obtained using density function (DFT) calculation [48].Our set of potentials reproduces lattice parameter to within 0.1%, that is, lesser than the other theoretical works.Available elastic properties of HAP materials are used to validate and refine the potential parameters.Table 3 shows the agreement between the calculated and experimental values of the elastic constants for the HAP structure.It can be seen that our obtained values are closer to those of experimental works and other previous theoretical works.The variability in the experimental elastic constants reported in the literature is shown in Table 3.
To determine Frenkel and Schottky type defect formation energies, isolated point defect (vacancy and interstitial) energies and relevant lattice energies were first calculated.Several possible positions were tested to confirm the optimal position of the interstitial site for defect occupancy, and the positions which have the lowest energy were taken for interstitial cations and oxygen, respectively.We also calculated antisite pair defects, which involve the exchange of a cation with another cation of different species.In addition, two other intrinsic defect schemes are considered.Scheme (i) involves one calcium vacancy and two hydroxyl vacancies, and scheme (ii) involves five calcium vacancies and two phosphorus interstitials.These defects are expressed by Kr ö ger-Vink notation and are shown in Table 4.The calculated defect formation energies, also listed in Table 4, were obtained by combining the energies of these point defects.The results indicate that the formation of defects involving phosphorus is generally unfavorable due to the high energy required to create this type of defect.The defects involving P ions are highly unfavorable due to strong P-O Coulombic interaction; that is, they would hardly occur compared to the other types of defects.The result also reveals a low energy (−0.71 eV) for scheme (i).This is consistent with the theoretical observation reported by de Leeuw et al. [31] using DFT calculation.The work of de Leeuw et al. [31] reveals that calcium vacancies compensated by hydroxyl are preferred to those compensated by phosphate vacancies.They also show that the calcium vacancies compensated by substitutional defects, such as carbonate groups, are more favorable for compensation by hydroxyl and that if the carbonate defect is accompanied by the substitution of a monovalent Na + or K + ion for the calcium ion, the defect formation energies are lower.In the more favorable schemes, calcium vacancies are involved, that is, one defect related to the calcium deficiency in the hydroxyapatite materials.Experimental works reveal that apatites are often found to be deficient in calcium compared to the stoichiometric material [51][52][53].Other mechanisms of defects, not considered by de Leeuw et al. [31], could be involved in the calcium deficiency in the hydroxyapatite.From our results, it can be seen that calcium pseudo-Schottky and calcium Frenkel defects are the next most favorable defects.In both schemes, calcium vacancies are involved.Thus, these defects also contribute to the higher concentration of calcium vacancy in HAP.This suggests that a higher concentration of calcium and hydroxyl vacancies could be present, which may contribute to calcium deficiency in the hydroxyapatite materials.
In the doping process, a description of the favorable substitution site and the charge compensation mechanism is very important information.From atomistic simulation it is possible to obtain quantitative estimates of the relative energies of different modes of dopant substitution.For this, trivalent (M 3+ ) and tetravalent (M 4+ ) dopants ions substituted at both Ca and P sites were investigated in order to help explain experimental results.The charge compensating mechanisms involved in the incorporation of M 3+ and M 4+ ions have not been clearly established, so all possible schemes are considered in this work.
M 3+ and M 4+ dopant ions can, in principle, be incorporated into the lattice at either Ca 2+ or P 5+ sites.In both cases, there is more than one possible mode of charge compensation.All possible charge compensation mechanisms were considered and are shown as follows.
Tetravalent dopants 2,5M 2 O 4 + 5Ca Ca + 2P P → 5M In the crystal structure of HAP [12], there are two types of Ca 2+ host site.Therefore, the incorporation in each site is considered.Seven reaction schemes have been considered for the incorporation of the trivalent ions and seven for the tetravalent ions.
In ( 1) to ( 4) the trivalent cation metals are considered at the Ca 2+ site, and in ( 5) to (7) the trivalent cation metals are considered at the P 5+ site.On the other hand, in (8) to (11) the tetravalent cation metals are considered at the Ca 2+ site, and in (12) to (14) the tetravalent cation metals are considered at the P 5+ site.In the all case, charge compensation is needed.
The solution energies for the trivalent defect were calculated by combining the appropriate defect and lattice energy terms and are listed in Tables 5-7.It is seen that the solution energy is different for the incorporation of the earth rare/metal transition with different ionic radius.The latter observation is explained by the lattice expansion/contraction, resulting in more/less space to accommodate earth rare/metal transition ions together with any charge compensating defects.From the tables it can be seen that substitution at the Ca site is more favorable for all trivalent dopants.The Ca site preference has been confirmed experimentally by Martin et al. [54] and Wakamura et al. [55].This behavior can be explained in terms of the difference in ionic radii between the trivalent dopant ions and the host site ions.The ionic radii of transition metal and rare earth ions vary around 0.55-0.75 and 0.86-1.03Å, respectively, and the ionic radii of Ca 2+ and P 5+ are 1.00 and 0.38 Å, respectively [56].The small difference between the trivalent dopant ions and Ca 2+ ions contributes to a small deformation in the lattice and consequently a small solution energy.Thus, the results show a degree of correlation between the dopant size and the solution energy.Dopants with a large difference in ion size from the host tend to be more energetically unfavorable.Between the two different Ca sites, it seems that, except for Lu, Er, Gd, Eu, Sm, and Nd show a preference for the Ca1 site substitution.
In the Eu-doped HAP, the 5 D 0 -7 F 0 transition of Eu emission is particularly informative.This is because the 5 D 0 state is unsplit, and the 7 F 0 ground state is also unsplit, so that transitions to it give straightforward information about the excited state.If more than one component is seen for this transition, it shows that there is more than one europium site [34].The two peaks attributed to 5 D 0 -7 F 0 transition were observed in the Eu emission spectra by Jagannathan and Kottaisamy [57] and Martin et al. [54], indicating that the Eu ion can be substituted at two sites.From our results, it can be noted that the solution energy for more favorable schemes for Eu substitution at the Ca1 site (1.79 eV) is close to that for the substitution at the Ca2 site (1.50 eV).This small difference justifies the possibility of Eu substitution at both host sites.The results of Martin et al. [54] support our hypothesis that Eu dopant ions prefer to be substituted at the Ca2 site due to the lower solution energy.They also report that the site occupancy ratio between the Ca2 and Ca1 sites is about 80%, with the majority of the Eu 3+ ions being in the Ca2 sites.
The charge compensation tends to occur basically via three schemes.In general, charge compensation via the formation of calcium vacancies is more favorable for charge compensation in the trivalent dopants, except for Er, Gd, and Eu, where the charge compensation is more favorable through           [58] reports that this value is below 1.67 in most cases and shows a decrease with increases in the La content in the samples; that is, it indicates the formation of nonstoichiometric HAP.This reduction of the Ca/P ratio is attributed to the formation of the calcium vacancy generated as a charge compensation defect needed in the substitution of the La ion in the Ca site.The occurrence of this scheme (see (2)) is the most favorable, as can be seen in Table 7. Similar behavior can be expected for all the trivalent dopants that are compensated by calcium vacancies.
In Table 8, the solution energies for the tetravalent defect calculated by combining the appropriate defect and lattice energy terms are listed.It can be seen that substitution at the P 5+ site by a compensated oxygen vacancy is more favorable for Mn 4+ and Cr 4+ dopants.This behavior also can be explained in terms of similar ionic radii between the tetravalent dopant ions and the host site ions.Mn and Cr dopants can exist as either M 3+ or M 4+ .It is necessary to know their precise oxidation state.According to our results, the favorable doping process for these polyvalent dopants should be an M 3+ form via the formation of calcium vacancies.

Conclusion
The basic defect chemistry including intrinsic defects and trivalent and tetravalent extrinsic dopants has been investigated using atomistic simulation.The most favorable intrinsic defect schemes are formed by one calcium vacancy and two hydroxyl vacancies.Trivalent and tetravalent dopant ion substitutions are found to take place preferentially on the Ca and P sites, respectively.In both cases, calcium vacancy defect formation is the more likely charge compensation mechanism for most of the trivalent and tetravalent dopants.In some trivalent ions, charge compensation by interstitial hydroxyl and interstitial oxygen is more favorable.The favorable doping process for these polyvalent dopants should be an M 3+ form via the formation of calcium vacancies.

Table 2 :
Lattice parameters (, , , , , and ) calculated in the present work compared to experimental and other theorical works.

Table 3 :
Calculated and experimental elastic constants of hydroxyapatite.

Table 4 :
Solution energy of intrinsic disorder in hydroxyapatite (eV/defect).

Table 5 :
Solution energies (eV) for substitution of the transition metals (M 3+ = Fe, Mn, Cr, and Sc) and rare earth (M 3+ = Lu and Yb) ions in HAP.

Table 7 :
Solution energies (eV) for substitution of the rare earth (M 3+ = Eu, Sm, Nd, Pr, Ce, and La) ions in HAP.

Table 8 :
Solution energies (eV) for different dopants in HAP.
[14]charge compensation for Eu ions is in agreement with that proposed by Martin et al.[54]and Ternane et al.[14].The production of stoichiometric HAP is obtained with a Ca/P ratio of 1.67.Experimental work realized by Mayer et al.