Simulation and Computer Study of Structures and Physical Properties of Hydroxyapatite with Various Defects

Simulation and computer studies of the structural and physical properties of hydroxyapatite (HAP) with different defects are presented in this review. HAP is a well-known material that is actively used in various fields of medicine, nanotechnology, and photocatalytic processes. However, all HAP samples have various defects and are still insufficiently studied. First of all, oxygen and OH group vacancies are important defects in HAP, which significantly affect its properties. The properties of HAP are also influenced by various substitutions of atoms in the HAP crystal lattice. The results of calculations by modern density functional theory methods of HAP structures with these different defects, primarily with oxygen and hydroxyl vacancies are analyzed in this review. The results obtained show that during the structural optimization of HAP with various defects, both the parameters of the crystallographic cells of the HAP change and the entire band structure of the HAP changes (changes in the band gap). This affects the electronic, optical, and elastic properties of HAP. The review considers the results of modeling and calculation of HAP containing various defects, the applied calculation methods, and the features of the effect of these defects on the properties of HAP, which is important for many practical applications.


Introduction
Hydroxyapatite (HAP) is a mineral compound from the apatite group, the basic salt of calcium and phosphoric acid (calcium phosphate) with the formula Ca 5 (OH) (PO 4 ) 3 , a widely employed multifunctional material in biomedicine, health care, biology, ecology, catalysis, and for environmental remediation [1][2][3][4][5][6][7][8][9][10]. First of all, HAP is the main mineral component of mammalian hard tissues (bone and teeth). Along with the organic component (collagen) and living bone cells (osteoclasts, osteoblasts, and osteocytes), HAP crystallizes within the gaps of stacked tropocollagen fibrils, forming and strengthening the bone structure [1][2][3][4]9]. Due to innate bioactivity and biocompatibility, HAP is a widely used material in medicine for bone and dental surgery. Here, hydroxyapatite is used as a filler to replace parts of lost bone (in traumatology, orthopedics, hand surgery) and as a coating for implants to promote new bone growth.
Since the mechanical properties of HAP differ from bone material, implants are made from the most suitable material (usually titanium is used), but its surface is covered with HAP [3], which provides better adhesion and reproduction of bone cells on its surface [2,4]. It is important to note here that biological HAP (bio-HAP) found in living organisms differs from an ideal synthetic mineral crystal of HAP in stoichiometric imbalance and various atomic PO 4 group [8]. Their formation is possible at the higher temperatures (of the order of~1100-1300 • C or about~1373-1573 K) [19,20] or also under radiation exposure [23].
Basically, the main physical and chemical properties of HAP are determined precisely by the presence of various structural defects, such as oxygen and hydroxyl vacancies, interstitials and substitutions of ions, atoms in the structure of HAP [1][2][3][4]. Various modern experimental and theoretical methods are used to study them. Among them, computer modeling and calculations, including first principles ab initio and density functional theory (DFT) methods [5][6][7][8][9] are the most important.
In previous research, using some modern DFT methods (including several types of the hybrid functionals and the theory of many-particle perturbations) and developed software tools (AIMPRO, VASP, and Quantum Espresso (QE) [7,8,[38][39][40][41][42][43][44][45][46][47][48][49][50][51][52]) in the modeling and computer studies of the properties of HAP structures with these various types of the defects, it was found that the forbidden band gap Eg of an ideal (defect-free) stoichiometric HAP lattice turns out to be much larger than the data of measurements of the optical absorption characteristics of the usual HAP samples [2,4,18,19]. These calculations show, that Eg can reach values of the order of~7.4-7.9 eV [7,35]. This significantly exceeds the experimentally observed values for the HAP samples used in practice (Eg~3.8-4.5 eV) [2,4,9,20]. These recent computational studies of the HAP structure with defects have convincingly shown that the deviations of the calculated band gap Eg of the ideal stoichiometric HAP crystal from the experimentally observed Eg values are caused precisely by the presence of the various defects in the HAP structure [2,[4][5][6][7][8][9]36,37].
The calculations of oxygen vacancies (from PO 4 and OH groups with different symmetry) and total OH-group vacancies in the HAP structure were carried out recently, showing very interesting results [5][6][7][8][9].
These calculations and studies have shown that, under equilibrium conditions, oxygen vacancies arise not only in the form of a simple vacant oxygen site (in a neutral charge state) [8], but can also be in the form of extended structures occupying several crystal fragments or chemical units (especially in a double plus charge state) [8,36,37]. In this case, the connection between the presence of defects and changes in the optical properties of HAP and the level of its photoelectronic work function is clearly traced [2,4,5]. In fact, here was shown, that for any different approaches and methods of DFT calculations, it is oxygen vacancies that play their decisive role here, changing the optical absorption width Eg of an ideal stoichiometic HAP crystal from its large characteristic values to a lower level of the order of~3.6-4.2 eV, determined by these introducing HAP defects (oxygen vacancies), which is as usually recorded on HAP experimental samples.
In general, these studies are still far from complete and require further deep continuation.
In this review, we present some recent data on the study of defects, primarily oxygen vacancies (of the various types) and the OH vacancies, obtained by modern methods of computer simulation and first principle of quantum-mechanical calculations based on DFT.
In addition, in this work, the calculations are compared with the available experimental data, and an analysis of this comparison is given.
One of the common structural peculiarities of HAP structure is connected with the pseudo-one-dimensional character of the apatite structure. Ca and OHions form a long chain along the main structural c axis [9][10][11][12][13][14][15][16][17][18][19], often being named as an 'OH-channel' in HAP. For the P6 3 /m phase consisting from two chemical formula units and the HAP general formula is [Ca 10 (PO 4 ) 6 (OH) 2 ], where the hydroxyl units OH shows stochastic orientation along OH channels.
Thus, according to X-ray data, this actually makes the material mirror-symmetric along the main axis, which runs along the OH channel. Conversely, when all OH blocks show the same alignment along the hexagonal axis, the mirror plane is lost and the space group symmetry drops to P6 3 . The ordering of the dipoles in the OH-channels, which interact with the chains of the OH-channel, also affects the property of the monoclinic phase. If OHions are oriented in parallel, HAP has non-centrosymmetric ordered structure (space group P2 1 ), which could reveal piezoelectric properties [17,18,36,37].
The structure of HAP, used to study the effect of defects on its properties, is primarily based on its initial pristine stoichiometric structural phase-hexagonal P6 3 , with a unit cell consisting of 44 atoms and containing structural OH channels with two hydroxyl OH groups in each unit cell [5][6][7][8][9][13][14][15][16][17][18][19][20]. Depending on the orientation of these OH groups, the cells can have different symmetry groups (Table 1): P6 3 /m-for the hexagonal disordered phase (when the orientation of the OH groups is random) and P6 3 -for the hexagonal ordered phase (when the orientation of the OH groups is parallel and directed in the same direction, which creates its own internal polarization, similar to ferroelectrics (see below, for example, in [36,37]). Some earlier calculations were performed using one unit cell model of hexagonal HAP P6 3 [5,6] (see Figure 1). The main peculiarity of the next study is the introduction of the supercells model made up of 2 × 2 × 2 = 8 HAP unit cells (space group P6 3 ) for hexagonal HAP phase both for pristine and defective HAP with OH and oxygen vacancies [7,8] (see below in Section 2.4).  6 (OH) 2 )-in P6 3 hexagonal phase. All OH groups are oriented in the same direction. They are positioned at the four corners of the unit cell, but only one pair in one corner belongs to this unit cell, the other three pairs belonging to neighboring unit cells (e.g., one OH per unit cell). It also shows the axes of the Cartesian coordinates x, y, z and unit cell parameters a, b, c. (Adapted with permission from ref. [5]; IOP Publishing, 2015).

HAP One Unit Cell Model
The stoichiometric HAP unit cell model for pure (defect-free) lattice in P6 3 hexagonal phase is presented on Figure 1 [5]. This model was used for the calculations firstly for LDA method [5,6] and then for GGA approach too in several cases [7]. It was used as an initial unit cell model of hexagonal HAP (P6 3 ) lattice and for the cases with several of its defects, such as, OH-vacancy, and various O-vacancies.
This one HAP unit cell consists of 44 atoms and contains two OH groups in each periodical unit cell. It should be noted that only one OH-channel with two OH groups is included in the unit cell, while the remaining three OH-channels (with six OH groups) belong to other periodically repeating unit cells.
For monoclinic HAP nanostructures, the unit cell consists of 88 atoms, because the lattice unit cell is doubled in the monoclinic case along the b axis. Besides this basic hexagonal unit cell of 44 atoms, the supercells totaling 176 atoms, made with two and four unit cells along the a and c directions of the hexagonal lattice were constructed to investigate defects in the work [5].
For the unit cell with 44 atoms, we used a 6 × 6 × 8 grid of k points (288 points), following the recipe proposed by Mokhorst and Pack [57], while for the large supercells, the grid used was 2 × 2 × 4. All k-point sets were folded according to the symmetry of the problem. All atom positions were allowed to relax along the forces acting upon them, until the total energy converged below 10 −5 eV. After convergence test calculations, a 2 × 2 × 4 k-point mesh was found good enough to sample the Brillouin zones (BZs) of both hexagonal and monoclinic cells.
The HAP unit cell model, containing several defects in the atomic positions of all the vacancies used in this work, is shown in Figure 2. Here, we indicate OH-vacancy and various types of O-vacancies in the HAP unit cell, as deviations from the initial stoichiometric hexagonal HAP structure [6]: (1) an O-vacancy in the OH group of an OH-channel structure (blue circle in Figure 2); (2) a full OH-vacancy from this OH-channel structure (purple circle in  For GGA approach, the Brillouin zone (BZ) was sampled using a Γ-centered 2 × 2 × 3 mesh of k-points. The Hartree-Fock exact exchange was evaluated at the same k-point grid used for the DFT potential and stored on a real-space grid of 128, 128, 96 points along a1, a2, and a3 lattice vectors, corresponding a, b, c (see Figure 1). The experimental lattice constants are a1 = a2 = a = 9.417 Å and a3 = c = 6.875 Å [11], corresponding to a grid density of about 14 points/Å along all three directions.
In all these calculations-except for the elastic constant calculations (bulk modulus B) and the DOS calculations-full relaxation of the cell, including changes in volume and adjusting the atomic positions and cell shape (total optimization of structure) was allowed. When calculating the elastic constants, only the relaxation of ionic positions within the strained cells was allowed. In several calculations, e.g., when defects as substitutions of atoms (Ca on Sr, Fe, Mg etc.) were explored, the calculations for fixed former atomic positions were also used, and then their relaxation to the new ones. The optimized data obtained for each HAp structure were taken (files with optimized and fixed coordinates for all unit cell atoms in atomic and Cartesian units were used) and transformed further to another, special HyperChem (*.hin) format (using the Babel program [56]). Then, all these atomic coordinate data files were uploaded into HyperChem software tools for further molecular modeling (with visual options for all studied HAp structures), calculations, and exploration of its physical properties using HyperChem package [40].

HAP Supercell Model
The main feature of our more detailed study is the introduction of the super-cells model made up of 2 × 2 × 2 = 8 polar HAP unit cells (with space group P6 3 ) for hexagonal HAP ordered phase (both for initial pristine stoichiometric and defective HAP with oxygen vacancies and with full OH group vacancy) [7,8].
It should be noted that this phase demonstrates polar ferroelectric state with the total polarization along c-axis, due to OH groups oriented parallel and along the OH-channels. While the disordered phase of hexagonal HAP (with space group P6 3 /m) forms non-polar paraelectric state, due to compensation of an opposite anti-parallel orientation of OH groups, which are along these OH channels [36,37].
Oxygen vacancies were introduced in hexagonal supercells made up of 2 × 2 × 2 = 8 polar HAP unit cells (space group P6 3 ), comprising a total of 352 atoms ( Figure 3). Although the monoclinic phase P2 1 /b was found to be more stable [12,33], they only differ by a few tens of meV per unit cell, with their electronic structures being essentially identical [34]. For the HAP unit cell, we have shown previously that convergence of the electron density and energy is obtained when sampling the band structure with a Γ-centered 2 × 2 × 3 mesh of k-points within the first Brillouin zone. By doubling the size of the cell along all principal directions, reciprocal lattice vectors are contracted by a factor of two so that a 1 × 1 × 2 k-point grid would actually improve on the sampling quality. From convergence tests, we found that the total energy of the 352-atom bulk supercell obtained with a 1 × 1 × 1 grid (Γ-point sampling) differs by less than 0.1 eV from a 1 × 1 × 2-sampled calculation. More importantly, relative energies and ionization energies of defects differ by about 1 meV only. Therefore, all defect calculations employed a Γ-point sampling [7,8].
A full relaxation of such large supercells using plane wave hybrid-DFT is prohibitively expensive. Instead, in works [7,8], defect structures were first found by relaxing all atomic coordinates within PBE, until the maximum force became less than 10 meV/Å. The resulting structures were employed on a second step, where the total energy was obtained within hybrid-DFT by means of a single-point calculation. This procedure was necessary due to the sheer size of the Hamiltonian at hand combined with the use of a plane-wave method.
As it was recently shown in [7,52], the relative energies obtained within this methodology are usually affected by error bars of the order of 10 meV or lower.

Main Structural and Mechanical Properties of HAP
We start by analyzing the structural and mechanical properties of HAP as initial stoichiometric as well as HAP with considering defects: OH-vacancy and various O-vacancy types. Tables 2 and 3 compare the calculated structural data (unit cell lattice parameters and volume) and bulk modulus (for some cases) with the respective experimental data. Table 2 presents the data of pristine stoichiometric HAP (computed by various methods and from different experimental data) and data of HAP with one full OH-vacancy (per one HAP unit cell). Table 3 consists of calculated data for various types of O-vacancies and it is shown below.   Figure 2). (2) OH-OH is the ordered P6 3 phase and OH-HO is the disordered P6 3 /m phase-comment 1 in the text.
The results obtained with LDA method (in earlier calculations using AIMPRO [5,6]) and with GGA approximation on the various DFT exchange-correlation functionals: PBE, HSE, B3LYP, and PBE0 (using VASP calculations [7,8]) are presented in Tables 2 and 3. The calculated lattice parameters for initial pristine stoichiometric HAP lattice within the PBE level approach are in good agreement with the various previously known PBE results (see [48][49][50][51] and references therein).

Comment 1 for OH groups orientation difference.
Recently computed data have shown, first, that HAP can coexist at different stages, both hexagonal and monoclinic [5,36,37], as well as in ordered and disordered hexagonal due to the orientation of the OH groups in one unit cell by aligning the OH units in opposite directions: OH-OH is the ordered P6 3 phase and OH-HO is the disordered P6 3 /m phase. The energy difference between the two phases is DE = E (P6 3 /m)-E (P6 3 ) [5,7,8]. This difference, obtained from our various calculations and presented here in Table 2, is: ∆E =~0.43 eV for PBE (GGA) and ∆E =~0.132 eV for LDA (AIMPRO) [5]. In a similar PBE calculation in [34], this ∆E~0.4 eV, and in the most developed model of a large supercell [7,8] is ∆E = 0.39 eV, and the last two energies were obtained from calculations of the completely relaxed B3LYP level [23].
Consequently, during heating (temperature above~1000 K) and cooling, many OH groups can change their orientation in different unit cells, which ultimately can lead to some stochastic redistribution of them, and the sample will contain a mixture of both types P6 3 and P6 3 /m hexagonal phase of HAP. It is very difficult to distinguish them experimentally. Therefore, the comparison of the results of our calculations and experimental data must be carried out very carefully, taking into account the possible errors arising also due to this factor of different random orientations of OH groups in experimental samples.

Comment 2 for two-steps calculation for supercell models.
For supercell models, structural optimization was made using PBE(GGA) only for initial stoichiometric HAP and after that the calculation in B3LYP was done in the fixed atomic positions in single point (SP) mode. All further calculations of the HAP with defects were performed using PBE(GGA) optimization only for atomic/ionic positions keeping the lattice parameters fixed as it was founded (established) for initial stoichiometric HAP lattice.
The results obtained show the usual~1% overestimation in relation to the experimental data. This is known to be mostly due to an artificial over-delocalization of the electronic density when the GGA approach is employed. It is shown that the calculations using B3LYP results in improvement of c lattice parameter, but parameter a is still overestimated by~1% [7]. This confirms previous reports that B3LYP generally overestimates the experimental lattice parameters as well [35]. On the other hand, our result differs from the previous B3LYP calculations of HAP, where a was underestimated by~1% [35]. More discussion about the reasons of these discrepancy was given in [7]. Table 2 also indicates that the lattice parameters calculated within hybrid DFT are generally closer to the experiments than those obtained using the semi-local functional. The best results are obtained for PBE0 with a deviation of <0.6% in relation to the experiments [7].
It should be noted, that the results obtained earlier in the LDA approximation are the closest to the experimental data on the lattice parameter a (even slightly better than in the case of a PBE0 with deviation less than 0.6%), while the value of the parameter c turned out to be more overestimated here (the deviation is almost 1.8%).
The bulk modulus (B) were obtained here usually by fitting the Birch-Murnaghan equation of state [5,[61][62][63][64]. The calculated bulk modulus show reasonable agreement with the experiments (see data in 2nd column of Table 2, Refs. [61,62]). The errors (shown in the table) were obtained from the standard fitting procedure. B3LYP calculations show deviation from the measurements by about 3% only, while other methods underestimate the experimental value of B by 7-8%. This level of accuracy is in line with typical discrepancies found for many other insulating materials [63]. It is interesting again to note, that LDA approximation shows here the same value of B as on PBE (GGA) level.
Experimental data on stoichiometric natural and synthetic samples should be used and analyzed rather carefully, since there may be impurity ions (carbonate, nitrate) that affect the lattice parameters. Therefore, the composition of stoichiometric hydroxyapatite should be confirmed by other analytical methods (for example, IR data), in addition to diffraction methods, which is usually used to determine the lattice parameters. Also, synthetic stoichiometric HAP should be well annealed at high temperatures in humid air. Table 2 provides experimental data that mostly meet these requirements [11,[58][59][60][61][62]. Below in Section 3.3.1, we discuss this data in comparison with our calculations.
The results of more detailed calculations for the presence of the various types of oxygen vacancies in HAP are presented in Table 3. Here, a complete optimization of the HAP structure (cell parameters and ion relaxation) was carried out only by methods LDA (AIMPRO) and PBE (GGA) in the one unit cell model, since such an optimization has not yet been performed for the supercell model.

Defects (OH-Vacancy and Various O-Vacancies) Influence on Structural and Mechanical HAP Properties
Let us now consider the effect of defects (oxygen vacancies and OH group vacancy) in the HAP structure, which changed in their structural and mechanical properties, and which turns out to be important for many practical and technological applications.
(1) Influence of OH-vacancy. The data calculated by LDA (AIMPRO) show that the relative changes in the lattice parameters of HAP after the appearance of an OH vacancy (in one HAP unit cell) both increase by small values: δa~0.15% and δc~0.05% (see Table 2).
At the same time, calculations using PBE (GGA) show that the value of the lattice parameter also relatively increases by δa~0.6%, and the parameter c increases by δc~0.5% (see Table 2). For comparison with experimental data, we use here the results of experiments on heating and cooling HAP samples, which should lead to the escape of some OH groups from the OH HAP channel and the formation of OH vacancies in HAP ( [58][59][60] in Table 2).
In situ diffractometric studies carried out in the group of Dr. N.V. Bulina showed that the cooling of apatite sample in helium (He) atmosphere, heated to a temperature of 1100 • C, leads to a decrease in parameter a and an increase in parameter c. Table 2 shows these obtained experimental data. In such conditions, OH vacancies formed at high temperatures remain up to room temperature and the cooled sample is now the oxyhydroxyapatite (OHAP) containing vacancies of OH groups. The changes obtained in the lattice parameters in the presence of OH vacancies in comparison with those known for pristine stoichiometric HAP turned out to be as follows: the relative changes were the largest for the parameter a, the decrease was up to δa~−0.016%, and for the parameter c, an increase by about δc~0.5%.
If we compare with the initial experimental data [11], then the changes are similar. Thus, these obtained data are very close to the calculated data in the main trend of changes in the lattice parameters during the formation of an OH vacancy: a slight decrease in the lattice parameter a (and in the calculations there is an increase) and a noticeable increase in c value for all cases due to the formation of an OH vacancy in the HAP sample. As for the volumes of the cells, with both methods of calculation, a relative increase in volume was obtained (Table 2). For Bulina's experimental data, the volume slightly decreased. However, if we compare it with the initial volume of [11], then it also increased.
It is important that the calculation results show a similar character of the cell parameter measurement. This means that the chosen model and calculation method are correct and allows us to investigate in more detail the physical reasons for this phenomenon [19].
The calculated values of the modulus B (performed after the LDA method optimized calculation) showed its decrease in this case of the presence of the OH vacancy in HAP sample, which means a decrease in the mechanical strength of HAP samples ( Table 2).
It seems that the presence of any vacancy type, including the OH group, should increase the lattice parameters, which is what we get. However, it is likely that in the experiment, out of two vacancies of OH groups, one disappears, and it is occupied by the O 2ion (that is, the hydrogen atom H also disappears here and one oxygen ion remainsthis is already forming now the oxyapatite (OAP or OA)). This is possible and leads to the following observed changes in the lattice parameters: HAP (a = 9.4236 Å, c = 6.8802 Å) goes into OHAP (a = 9.4155 Å, c = 6.8835 Å). Moreover, some recent experimental data show that if the OAP is formed here, then the lattice parameters should be as follows: OAP (a = 9.4057 Å, c = 6.8938 Å).
Unfortunately, we have not yet performed calculations of the exact OAP models, containing two OH vacancies per one HAP unit cell. It would be very interesting and useful, and we will definitely do it in the near future. In this case, we suppose to reach decreases of both lattice parameters as experimentally observed by Bulina.
(2) Influence of only one O-vacancy from OH group. In this case (see in Table 3), the calculated values of the HAP unit cell parameters give some relative decrease in the parameter a (both for the LDA (AIMPRO) and GGA (PBE/VASP) methods), while the c parameter also decreases when calculated by the LDA method, but slightly increases when calculated by the GGA (PBE) method. Unfortunately, we have no experimental data for this case.
(3) Influence of the different O-vacancy from PO4 group. In this case (presented in Table 3), the calculated values of the unit cell parameters were considered for different positions of oxygen atoms that create O-vacancies in various place in the HAP unit cell (for example, see Figure 2). For the model of one unit cell, these positions were selected and calculated by the LDA method in the works [5,6], and in Table 3 these data are given. Similarly, the calculation using GGA (PBE) method was performed and presented here too.
These O-vacancies were computationally investigated in more detail by the GGA (PBE) and B3LYP methods in the work [8] in the HAP 2 × 2 × 2 supercell model. We will analyze these data separately in Section 3.3 below, devoted to supercell model (especially devoted to its electronic and optical properties).
For the cases of O-vacancies from the PO4 group indicated in Table 3, the results obtained by both LDA and GGA (PBE) methods in comparison with initial stoichiometric HAP (hexagonal P6 3 ) from Table 2 show that both the unit cell parameters a and c decrease (as well the unit cell volume), despite some of their differences for the different positions of these selected oxygen atoms, create these O-vacancies. We do not have experimental data for structural analysis in this case, but these cases demonstrate very interesting and important results of their electronic and optical properties, which will be analyzed below.

Electronic Properties of Pure HAP by Various Methods
To study the electronic and optical properties of HAP, both the initial defect-free and those containing various defects, the densities of states of electrons (DOS) at the different energy levels were calculated. This makes it possible to determine the effect of various defects on the energy band structure of these HAP properties.
As a result of LDA calculations (using AIMPRO software [38]) and GGA calculations (using VASP with functional PBE [39,48]), the distributions of the DOS corresponding to the filling of energy levels in the band structure with electrons and the energies of HAP according to the band theory were obtained [5][6][7][8][9]36,37]. This makes it possible to determine the main energy band structure's parameters-the position of the top of the valence band Ev, the bottom of the conduction band Ec, and the value of the forbidden band gap Eg = Ec − Ev, as well as the position of additional energy levels Ei induced by defects in the internal range of the Eg.
The question now is how defects change this entire band energy structure of the DOS and how close these changes it will be to the experimental values, when creating certain defects in the HAP structure. For these calculations at the first stage [5,6], the defects simulated were created by vacancies into one hexagonal unit cell of HAP consisting of 44 atoms (Figure 1). Figure 2 schematically shows examples of simulated HAP structures with defects such as oxygen vacancies from the OH and PO 4 groups, as well as a vacancy of the entire OH group. For more details and precise calculations with various DFT hybrid functional we use supercell model consisting of 352 atoms (eight unit cells) (Figure 3) [7,8].
Let us first consider HAP without any defect in the initial pristine stoichiometric hexagonal phase P6 3 symmetry. Figure 4 (Left figures: a,b,c) shows examples of DOS for perfect stoichiometric defect-free HAP (taking into account both deeper energy levels in the valence band and in the vicinity of the forbidden band) obtained in LDA calculations [5,6]. Similar GGA (PBE) calculations were carried out using VASP software [39]. Figure 4c shows results of DOS from GGA (with PBE functional). Figure 4a-c denotes Eg 0 = Eg because it refers to a defect-free HAP (see details also Section 3.2.3 below).
Note that the calculations using the LDA and GGA (PBE) approximations generally do not fundamentally differ here, the difference is only in the energies-GGA (PBE) gives a larger value of the band gap Eg~5.26 ± 0.05 eV (for perfect stoichiometric HAP) compared to calculations by the LDA method, which give band gap value Eg~4.6 ± 0.05 eV. These data are close to results of other authors [30][31][32][33][34][35].
Further, Figure 4 (Right figures: (a−e)) shows the total DOS of a HAP unit cell obtained within PBE to be more precise with detailed calculations [7,8], including also local density of states (LDOS) [7].  Here, we can see that the calculations using both GGA (PBE) approach do not much differ between one another, the difference is only small in the energies; first GGA (PBE) gives a value of the band gap Eg~5.26 ± 0.05 eV, while second PBE gives a value Eg = 5.23 eV, that is in the accuracy frame and in line with calculations of other authors by this method [9,[30][31][32][33][34][35][36][37].

Electronic Band Structure of Perfect Stoichiometric HAP
Further calculations of HAP perfect stoichiometric structure and properties were developed and performed in [7] using various DFT hybrid functionals. These results obtained demonstrate rise of Eg for perfect stoichiometric HAP using a more developed functional and method. We discuss this point below.
Interestingly, the energy scale of above DOS calculations is directly comparable to the band structure obtained in [7] calculations (see Figure 5a for PBE). Figure 5a-d compare the electronic band structure obtained using different exchangecorrelation functionals (PBE, HSE, B3LYP, and PBE0) with the analogous G0W0 quasiparticle calculation shown in Figure 5e [7]. The band energies along the several highsymmetry directions were obtained by interpolation of the first-principles data using Wannier90 [65].
The shape of the PBE band structure in Figure 3a is indistinguishable from that reported by Slepko and Demkov [34] displaying a low-dispersive valence band top and a high dispersive conduction band bottom (thick bands). Dispersion of the conduction band minimum states is considerably more pronounced along directions parallel to the c-axis (Γ-A, K-H, and M-L), indicating a stronger carrier delocalization and mobility along the main axis. This property could be explored for tuning HAP electrical conductivity through n-type doping or for photo-current measurements. On the other hand, p-type doping is not expected to be beneficial. The valence band top states show very little dispersion, and their heavy holes imply a relatively lower mobility.
Also, in agreement with [34] we find HAP to be an indirect-gap material with Eg = 5.23 eV at the PBE level. The conduction band minimum is located at k = Γ, while the valence band top energy was found somewhere along Γ-K or Γ-M. The valence band maximum along Γ-M is only 0.1 meV higher than the one along~Γ-K. We note that this picture was the same regardless of the functional used, including when using the G 0 W 0 method.
The band structure obtained within the HSE06 level is shown in Figure 5b. The increase in the band gap width by more than 30% with respect to the PBE result is selfevident. Using HSE06, we obtain Eg = 7.11 eV. Often, the band energies are offset in order to lock the valence band top at the origin of the energy scale. We did not follow this procedure and that allowed us to disclose how the gap change depends on the shift of both valence band and conduction band states. Figure 5b-d show that admixing a fraction of Fock exchange with the semi-local exchange energy has a profound effect on both valence and conduction band states. Consequently, the use of hybrid functionals has implications not only to the accuracy of calculated defect-related or inter-band transitions (e.g., observed in luminescence or UV-VIS absorption), but also to transitions involving core or vacuum states (e.g., observed in electron photoemission or core electron energy loss spectroscopy) [7].
The gap of the B3LYP band structure depicted in Figure 5c is 0.6 eV narrower than that reported by Corno et al. [35] using the same functional. Again, the difference is likely to be due to the unsuitability of the atomic-like basis employed in [35], which resulted in the under-screening of the band structure. The band gap width and the band gap edge energies obtained at the B3LYP-level are closer to the G0W0 results than any other functional.
Thus, here it should be noted and concluded that with the advent of each of the next developed types of DFT functional, the total calculated electronic band structure of pure HAP is change: the energies of valence band top is down and the bottom of conductive band is up. As result, the width of the forbidden band gap Eg = Ec − Ev increases (see Figure 5a-e).
This trend of the calculated band gaps agrees with that obtained by Garza and Scuseria [66] for an eclectic mix of semiconductors and insulators. Accordingly, HSE06 and B3LYP showed a closer correlation with the experiments (the former giving slightly smaller gaps overall), whereas inclusion of larger fractions of exact exchange, like in PBE0, led to an overestimation of Eg. This ordering also suggests that the gap width obtained within G 0 W 0 (Eg = 7.4 eV) should be close to the real figure.
Another interesting point should be noted here, regarding the origin of the bottom of the conduction band, in particular, the highly dispersive bands shown in Figure 5, the situation here is somewhat controversial. In [34], an analysis of LDOS at the bottom of the conduction showed that the lowest energy bands originate from Ca(4s) states.
The authors in [7] argue that this view finds several difficulties. First, in [7], the authors find the onset of the Ca-LDOS (E onset = 8.40 eV) located~1 eV above the conduction band minimum energy. Within that energy range, all that can be related to Ca are the flat Ca(3d) bands above 8.4 eV. Second, right Figure 4e shows that states just above Ec have a considerable localization on OH molecules.
Interestingly, in this case the highly dispersive bands that form the bottom of the conduction band of HAP are anti-bonding states from an infinite . . . OH-OH-. . . hydrogen bridge sequence (along OH-channel through all HAP crystal structure), much like a 1Dice phase. This statement finds support in the band structure of hexagonal ice (see, e.g., Figures 4 and 6 in [67,68]). From comparison of these band structure images with Figure 5a, it becomes immediately evident that the dispersion shape of the lowest conduction band of HAP is analogous to that of the hexagonal ice.
In this relation, it is interesting to note that nano-confined water [66] inside hydrophilic cavity (or inner channel) inside so-called peptide nanotube (such as, e.g., the diphenylalanine nanotube) has structure close to the one-line hexagonal ice structure [69,70]. So, we should note here that it is a common property in tubular or columnar molecular-like periodical crystal structure, which could be observed in these related structures.

Electronic Band Structure of Defective HAP (with OH-Vacancy and O-Vacancies)
In Figure 6, examples of DOS in the results of GGA (PBE) calculations are given for defects in HAP structure, such as O vacancies from the OH group, O vacancies from the PO 4 group, and a full vacancy of the OH group, based on HAP unit cell model with 44 atoms.
These data are close to results of other authors [19,[30][31][32][33][34]. The results obtained also lead to a small shift in all level energies for defects in the calculations of the GGA (PBE) methods relative to LDA (Table 3). Both methods predict changes in the band gap Eg and the electron work function ∆ф, which can be measured experimentally [2,4,9,53]. To compare this change in the band gap Eg in the presence of a defect in the HAP and to distinguish it from defect-free HAP, we here introduced the designation for the case of a defect-free HAP as Eg0, and then the relative change in the band width will be ∆Eg = Eg − Eg0, where Eg is the band gap width in the case of a HAP with a defect. The electron work function фhas approximately the same changes ∆ф≈ ∆Eg (it is important for comparison of such changes with experimental data of the photoelectron emissions) [4,5,9]. All these computed changes are summarized in Table 4 below. In addition, we will denote the energy of photoexcitation of an electron from the local level Ei to the conduction band Ec, which arises in the presence of a defect, as Eg* = Ec − Ei, and this will correspond to the new effective width of the band gap Eg* for optical excitation of an electron and/or photoabsorption of light, which are usually observed in the experiment. These designations are also used in Figure 6 and others figures and tables below.
Similarly, here we introduce the excitation energy of an electron from the valence band Ev to this level Ei, as Ev − Ei, which is presented in Table 4.
It is important to note that the levels of optical absorption for OH vacancies in HAP, established here (Table 3), lie in the region of~1.6-2.3 eV and they are practically not observed in the experiment, even at a significant concentration of OH vacancies. The calculated density of electronic states (DOS) from the HAP data shows that the intensity of the DOS peaks for the energy levels Ei caused by OH vacancies (located approximately in the middle of the band gap Eg) is very low, compared with the intensity of the main DOS peaks (the top of the valence band Ev and the bottom of the conduction band Ec, which determine the band gap Eg = Ec − Ev) [5]. Therefore, this transition is difficult to detect experimentally by photoelectron excitation. This turns out to be possible only by the photo-luminescence method, upon excitation of a significant number of electrons and from deep valence band levels into the conduction band by the method of synchrotron irradiation [53]. However, in this case, namely the changed value of the optical absorption can be experimentally recorded, which rise up to the value Eg*~5.3-5.5 eV, due to rise of HAP band gap caused by OH vacancies (Table 3).
At present, work on the experimental detection and observation of these spectral characteristics of optical absorption in this energy range is underway, created by such defects, as the vacancies of OH groups in HAP, are being actively carried out and are in progress. Figure 7, Figure 8, and Figure 9 below show examples of the band structure and the partial charge density for electrons with their eigenvalues in the range specified around the energy level Ei of corresponding defect structure for main defects in HAP, obtained after optimized VASP calculations in one unit cell model.   These results clearly show that it is precisely the oxygen vacancies of the PO 4 group that provide the optical band with the energy Eg* =~3.6-4.2 eV, which is usually observed in many experimentally prepared samples [4][5][6][19][20][21][22]36,37,53]. All these samples in the manufacturing process undergo heat treatment at temperatures T = 700-1200 • C (that is 973-1473 K) and this leads to the formation of not only OH vacancies, but also a multitude of O vacancies [19][20][21][22] (especially from the PO 4 group). This leads to the observed optical absorption and optical band width Eg*~3.6-4.2 eV, while the band gap of the defect-free perfect stoichiometric HAP turns out to be wider and equal to Eg~5.4 eV and more (in the GGA (PBE) calculations).
Note that the existing various experimental data [6,9,20] demonstrate optical absorption in HAP samples just in this energy range. In addition, these calculation data correlate with other, more accurate calculations using more advanced methods and models, considered below in the Section 3.3 (Section 3.3.1, Section 3.3.2). This means that even simpler models and methods correctly reflect the essence of changes in the physical properties of HAP in the presence of defects in its structure.
Further development of the calculations was the transition to calculations from one unit cell to a super-cell model consisting from of 2 × 2 × 2 = 8 unit cells, (space group P6 3 , comprising a total of 352 atoms), as well as the application of the hybrid functional PBE (for 1st step optimization) in combination with B3LYP (for2nd step calculation) using DFT method [7,8]. This made a possibility to classify different types of oxygen vacancies of the PO 4 group more correctly and accurately, which increases the calculation accuracy and to highlight not only neutral defects, but new more complex types of defects-extended charged complex oxygen vacancies.
These new structural extended (or bridging structures) defects can play an important role in changing the properties of electrically charged HAP.

Oxygen Originated Complex HAP Defects in Supercell Model
Performing these series of calculations, it was found that different types of defects arise, if the symmetry of the atomic group associated with different spatial arrangement of various atoms in the PO 4 group. Besides, usual oxygen vacancy V O , here arises complex extended defects, which depend on their charge state Q = 0, +1, +2. (Figure 10) [8,36,37]. The O vacancy of OH group influences the formation of new complex defects because of charge Q variation. Crystalline HAP has oxygen atoms on four symmetries in equivalent sites, which are referred to as oxygen types I−IV. Types I, II, and III in the phosphate units, while type IV oxygen atoms are located in the OHanions. These are denoted as O(I), . . . , O(IV) and are shown in the upper half of Figure 10, where portions of bulk HAP are depicted. We note that PO 4 groups have two nearly symmetric O(III) atoms which are superimposed in the upper view of bulk HAP in Figure 10.
Among the many vacancy structures investigated, those shown in Figure 10a−e are the most relevant as they showed lower energy. Metastable structures with more than 2 eV above the ground state (for each particular charge state) will not be discussed. Figure 10a describes the formation of a pyramidal PO 3 structure, where removal of the O(I) atom in "Bulk HAP" highlighted using a bright red color, leads to V O(I) as shown in the lower part of the figure, where a "Defective HAP" region is shown.
Analogous structures for V O(II) and V O(III) were obtained as well. In the neutral charge state, the resulting [PO 3 3− ] PO4 0 structures display a fully occupied sp 3 orbital on the P atom, resembling the phosphine molecule. This is shown in Figure 11a, where an isosurface of the electron density corresponding to the highest occupied state is represented in blue for the specific case of V O(III) 0 . Reprinted with permission from [7]. (Adapted with permission from [7]. CCC (2018) AIP Publishing). Figure 11b shows the case of a missing O(IV) atom, leaving an isolated H atom in the OH channel. This defect can be rationalized as the removal of neutral O from OH − , leaving a hydride anion, and maintaining charge neutrality of the whole system. In the neutral charge state, we have a [H − ] OH 0 structure, where after structural relaxation, the hydride species becomes located close to the site of the missing O(IV) atom (compare lower diagrams of bulk and defective HAP in Figure 10b).
The density corresponding to the highest occupied state of V O(IV) 0 is represented in Figure 11b, which clearly shows the formation of a hydride anion in the OH channel. Figure 11e shows the OH-vacancy defect [7]. In V O(I) −V O(IV) defects, all atoms (except the missing oxygen) remain essentially close to their original crystalline coordinates, hence the use of subscripted O(I)−O(IV) labels to identify their structure. However, additional structures, hereafter referred to as extended structures, were also found for the oxygen vacancy in HAP [8]. One type of such extended structures can be described as a pair of neighboring oxygen vacancies connected by an O-interstitial, 2V O + O. Another 2nd type is best described as a complex made of an OH-vacancy next to an H-interstitial, V OH + H. Two defects of type 2V O + O are singled out and labeled V O(A) and V O(C) . They are shown in Figure 10c,e, respectively. One V OH + H defect is shown in Figure 10d and is referred to as V O(B) . These new defect structures are very important for practices too.
A dashed circle is used in the figure to highlight the missing OH unit. The highest occupied state of the extended structures is shown in Figure 11c,d. They either overlap the void regions of the HAP crystal or the vacant volume of the OH-channel, thus suggesting that they are donors with antibonding character or strong resonance with conduction band states [7].
For relaxations in charge state +2 that started from structures I and III, the final structures were, respectively, A and C. Here, the P atom of the PO 3 unit in V O(I) (or V O(III) ) moved across the plane defined by the three O atoms to connect to the O atom from the nearest PO 4 moiety. Such a severe relaxation can be explained by electron transfer from a neighboring PO 4 3− anion to the empty P(sp 3 Figure 10d. The study of transitions between the energy levels of such more complex defects (including at different levels of the charge state of the HAP crystal, and different types of its conductivity: electron n-type, or hole p-type) showed that the transition of electrons can also include a phononless transition with energy~3.6-3.9 eV [8,36,37]. This mechanism can also explain the onset of absorption at 3.4-4.0 eV in the experimental observation of photocatalysis in HAP under constant UV illumination [6,8].
The results obtained convincingly show the primary role of defects such as oxygen vacancies and vacancy of OH group in the formation of actually observed values in the range approximately of~3.6-4.3 eV for oxygen vacancies Ovac {VO(I), VO(II), VO(III)} from PO 4 groups, and at~5.06 eV for oxygen vacancy VO(IV) from OH group (see below in Section 3.3.3) for the effective optical absorption and excitation bands Eg* in various samples of HAP materials [5][6][7][8][9]. At the same time, for example, the main experimentally observed values of the optical absorption band according to data measured UV spectra for the HAP samples of marine origin, and the commercial ones, respectively, (presented in Tables 2 and 4 in work [6]) lie in the close range approximately of~3.4-4.2 eV (which could be arisen by oxygen vacancy from PO 4 groups) and~5.05 eV (which could be arisen by oxygen vacancy from OH group). These data directly correspond to the data calculated in our computational studies presented and discussed in this article.
Nevertheless, the issues related to obtaining more accurate, fair, and correct numerical results, which are very important for many practical applications, as well as the issues of statistics (and dynamics) of accumulation and annealing of defects in various unit cells in the entire volume of the crystal, remain rather complex problems and require further research.
In addition, it is obvious that all these points require additional experimental observations and research. This work must be carried out in a complex.
These works are now in progress in the computational direction, and new results will soon be obtained in this direction of calculating defects in HAP, including, in addition to the types of vacancies here considered, also other defects of the type of insertions and substitutions of some atoms in HAP unit cell lattice.

Kohn−Sham Energy Levels of Neutral Oxygen Vacancy V O Defects
Inspection of the Kohn−Sham eigenvalues at k = Γ confirmed that neutral vacancies V O defects are all donors [8,37]. The calculations performed in these cases, taking into account the even more accurate B3LYP calculation scheme, showed similar shifts in the energy levels created by defects such as oxygen vacancies of the PO 4 group, and their main contribution, which determines the change in optical properties and the change in the work function, remains at a level of~1 eV or less (Table 5), despite the differences in the calculations of the band gap Eg for initially perfect stoichiometric defect-free HAP. As can be seen, regardless of the method of calculation, oxygen vacancies of the OH and PO 4 groups (in the absence of charge in these defects, i.e., at Q = 0) form a group of energy levels ( Figure 12) located close to the top of the valence band and they are donor electrons [8,48] (electron acceptors, that is, levels with a negative charge were not found). In this case, the shift of the levels of oxygen vacancies in the PO 4 group is~1.15-1.65 eV upward from Ev, and for a vacancy from the OH group,~0.4-0.7 eV, and this also corresponds to a change in the electron work function during the formation of such defects. As everyone can see, the positions of these energy levels are not very different (especially in the case of PBE calculations), although in the case of B3LYP these levels are slightly higher than Ev in the direct case. In any case, these deviations are within 1 eV.
It is also important to note that such levels close at the top of the valence band have recently been observed in experiments on photoelectron emission spectroscopy, and the work function of a photoelectron from HAP was measured for various external influences [2,4,5,53]. In addition, in [53] it was noted that such energy levels (and photoelectrons emitted from these energy levels) can arise under a number of actions on HAP samples (heating and annealing, gamma irradiation, microwave effects, and combined hydrogenation with microwave radiation).
In these cases, a sufficiently large number of oxygen vacancies (also in the OH group) can be induced having the lowest energies levels measured from the top of the valence band (see Table 4 and Figure 12).

OH-Vacancy in HAP Supercell Model and Some General Remarks
The vacancy of one OH group in the HAP supercell model was considered in a similar way. In this case, calculation with the GGA approximation in the hybrid functional PBE and B3LYP was performed using VASP [39] and Quantum ESPRESSO [41]. Figure 13 shows the results obtained in this case. After PBE optimization and B3LYP calculation using VASP energy band structure for HAP with OH-vacancy (Figure 13b), the following data are presented: Ev = 0.726 eV, Ec = 8.066 eV and Eg = Ec − Ev = 8.066 − 0.726 = 7.34 eV; the defective energy level due to OH-vacancy Ei = 5.3586 eV (occupied, spin up) and resulted optical gap from this energy level is equal Eg* = Ec − Ei = 8.066 − 5.3586 = 2.7074 eV. This data is in line with previous data obtained using another similar method in one unit cell model. In addition, it is important to note that the energy of formation of defects in the HAP crystal [7,8] was considered and investigated also in the formalism with a chemical potential based on work [52], which describes the procedure for calculating the energies of formation of defects in a similar study. This formalism also takes into account the formation of defects with an excess of electrons ∆n e (with respect to the neutral state) when they are capable of capturing/emitting electrons from/to an electronic reservoir with a chemical potential µ ε . The results obtained in this approach allow more correct calculations taking into account the thermodynamics of the formation of a defect with a screened charge as a quasi-particle in a crystal [7,8]. As a result, we obtain corrections to the energies, including the energies of optical properties and transitions, which, in our opinion, turn out to be even closer to those observed experimentally. These generalized data are shown here in Table 6. from OH} as well as V OH vacancy of OH group-in supercell model; for one unit cell model vacancies Ovac notation are the same as in Table 3 Ovac from Figure 2 Ovac from Figure 10 Kohn− Sham The results obtained convincingly show the primary role of defects such as oxygen vacancies and vacancy of OH group in the formation of actually observed values of the effective optical absorption and excitation bands Eg* in various samples of HAP materials. These works are in progress and new results will soon be obtained in this direction of calculating defects in HAP, including, in addition to the types of vacancies here considered, also other defects of the type of insertions and substitutions of some atoms in HAP unit cell lattice.

Conclusions
Computer studies of defects in hydroxyapatite, carried out by various methods, show their significant influence both on the structural and mechanical, as well as on the electronic and optical properties of hydroxyapatite with such defects. This is very important for practical applications. In this review, we consider in more detail defects such as OH group vacancies and oxygen vacancies from different PO4 groups and OH groups as well. Calculations performed by different methods showed that these vacancies significantly change, first of all, the electronic and optical properties of HAP. In this case, the vacancy of the entire OH group leads to an entire absorption band in the range~1.4-2.4 eV (with several DOS peaks). It is close to the red color of optical spectrum. However, the intensity of these DOS is small compared to the DOS of the top of the valence band. As a result, it is practically impossible to register this range of optical absorption (it is at the level of measurement errors). At the same time, an OH vacancy leads to a noticeable increase in the band gap of the order of ∆Eg~0.5-0.9 eV. This leads to a change in the optical absorption and also a change in the work function of the electron exit, which is recorded in the experiment. For example, it can be the level Eg~5.5-5.75 eV, which belongs to the ultraviolet spectrum and is close to the observed optical absorption of HAP samples with OH vacancies (created in them experimentally upon heating and subsequent cooling).
Oxygen vacancies have very different optical properties. Oxygen vacancies arising from various oxygen atoms of the PO 4 group and the OH group turn out to be of different types and depend on the symmetry of the corresponding group. As it was shown in [7] that besides neutral O-vacancy here may exist charged oxygen vacancies also form not only point defects, but more complex defects-extended or bridging (in the case of a charge Q = +2, Figure 10). The latter reconstructing to point-like defects (at Q = +1 and = 0). This rearrangement from extended to point defects occurs with bond breaking and causes also optical absorption effects. This transition leads to spontaneous rupture of bridging P-O-P or O-H-O bonds at extended defects and, most likely, explains the onset of absorption at 3.4-4.0 eV for observing photocatalysis under constant ultraviolet illumination [6]. It is important that these types of structure and the stability of defects strongly depend on the charge states.
It was found that the oxygen vacancies essentially occur in such two distinct forms, either as a simple vacant oxygen site (referred to as structures I-IV, Figure 10a), or as an oxygen atom replacing two neighboring oxygen vacancies (bridge or extended structures named as "A-C", Figure 10c-e). The former type of vacancies is deep donors, while the latter are shallow donors with rather low ionization energies. No acceptor states (stable negatively charged defects) were found. Vacancy structures I-IV are more stable in the neutral charge state, while bridge structures A-C are preferred in the double plus charge state. This means that the oxygen vacancy adopts rather different configurations on samples where the Fermi energy is in the upper or the lower half of the band gap. As regards to the neutral O-vacancy, corresponding inspection of the Kohn-Sham eigenvalues at k = Γ confirmed that neutral oxygen vacancy defects are all donors involving the luminescence and absorption of~3.6-4.2 eV.
Further development and more accurate calculation of these electronic properties and optical photoexciting and photocatalytic processes can be made also by addition-correct calculations of the electron-electron correlation of the excited electron states taking into account the Frank Condon (FC) relaxation. In [36,37], these contributions are calculated for the case of the V O (IV) vacancy in HAP P6 3 supercell model and their results showed that the FC shift obtained has reasonable values. These investigations should be continued for all other oxygen vacancy types.
Despite the need to continue further comprehensive study of HAP defects, however, based on the results already obtained so far, we must conclude that for any case, irrespective of these FC relaxation processes, the formation of the various types defects in HAP through different oxygen vacancies gives rise to an opportunity of the photo-excitation processes in the close ultraviolet (UV-C) and in the visible light region. These photo-excitation effects obviously produce the photo-catalytic activity of HAP, as well as provide the changes in the photoelectron work function and surface electric potential, which is important for implant covered by HAP. One only needs to introduce a sufficient amount of these oxygen vacancies by some external actions, such as heating/annealing, gamma-irradiation, and combined hydrogenation with microwave irradiations [71]. The electron photoexcitation from additional energy levels (with Eg* =~3.4-4.2 eV), which arose due to oxygen vacancies inside the forbidden band gap, can provide this electronic and optical properties of the treated HAP samples.