Using DFT to Calculate the Parameters of the Crystal Field in Mn2+ Doped Hydroxyapatite Crystals

Crystal field parameters for two nonequivalent positions Ca (I) and Ca (II) for hydroxyapatite (HAp) crystals from the density functional theory (DFT) are calculated. Calculations are compared with the experimental electron paramagnetic resonance (EPR) spectra (registered at two microwave frequencies) for the synthesized Mn-HAp powders Ca9.995Mn0.005(PO4)6(OH)2. It is found that in the investigated species, the manganese is redistributed between both calcium sites with prevalence in Ca (I). Agreement between the calculated and experimental data proves that crystal field parameters in HAp can be calculated in the classical DFT model using the distributed electron density.


Introduction
Hydroxyapatite (HAp)-based (bio)materials are widely employed for tissue engineering applications, especially in orthopedics and dentistry [1]. The crystal structure of HAp (Ca 10 (PO 4 ) 6 (OH) 2 ) allows for a wide range of substitutions, which can affect various material properties ranging from mechanical to antimicrobial [2][3][4]. Among various ions with therapeutic effects, manganese (Mn) has recently attracted attention for inclusion in diverse bioceramics [5][6][7][8][9]. Manganese is one of the essential trace elements for human health and one of the metals found in natural (bio)apatites. The presence of manganese in the HAp structure, for example, is able to change the adherence of bone cells to an implant material, favors osteoblast activity and proliferation, and enhances antimicrobial properties [10]. Manganese can induce alterations in mineralogical structure, physical properties, reactivity, and solubility of HAp, and influence the magnetic properties of a substance [11]. Hence, it is important to understand the incorporation mechanism of Mn into HAp, as well as structural alterations with Mn contents [11,12].
Although the incorporation of Mn into the HAp crystal structure has been extensively investigated, the preferred Ca site for Mn substitution is still under debate. There are two structurally distinct Ca sites for manganese localization in the HAp structure [13,14]. The Ca (I) site is nine oxygen coordinated, while the Ca (II) site is surrounded by six oxygen atoms and one hydroxyl group. It was shown, for example, that the localization and redistribution of Mn 2+ between the calcium positions in the inorganic component of the human aorta atherosclerotic plaques (composed mainly of the carbonated, nonstoichiometric hydroxyapatite) influence plaque stability and, therefore, can serve as a surgical indication [12,15].
Electron paramagnetic resonance (EPR) is known as one of the analytical, nondestructive methods for determining the type and concentration of paramagnetic impurities in liquids and solids including HAp [16]. To calculate the concentration, the integral signal intensity should be known [17]. To analyze the EPR spectra in more detail, the spin Hamiltonian is used, which, in general form, is the sum of the Zeeman interaction, zero-field splitting parameters (ZFSs), and hyperfine interaction. The fine interaction reflects the influence of the electric field created by other ions and electrons on the level splitting of a paramagnetic center with an electron spin S > 1/2. The corresponding spin Hamiltonian can be written in the formĤ ZFS = ∑ k=2.4.6 ∑ k q=−k B q kÔ q k , where B q k are the parameters of the electric crystal field, andÔ q k are the Stevens operators. Stevens operators are known, while the parameters of the crystal field depend on the environment of the impurity ion. It is possible to calculate these parameters for various positions of the impurity center, simulate the EPR spectra of these centers, and compare the calculations with the experimental data. Thus, it gives an ability to draw conclusions about what position is most likely for the impurity center. However, the EPR spectrum of powder is an average EPR spectrum of a large number of microcrystals differently oriented that prevents simple decoding of the center position. Therefore, in our work, we propose a calculation method that helps identify the positions of the impurity center.
To calculate the parameters of the crystal field, quantum mechanical approaches are most often used [18,19]. However, the classical approach can also be used, such as calculation of the crystal field parameters in the model point charges through the Coulomb interaction. In this work, a combined approach is used, where the electron density distribution is calculated using the density functional theory (DFT) method, and then the parameters of the crystal field are computed using the calculated electron density and atomic cores (at points of interest).
Nowadays, with an increase in computing power and a growing interest in materials based on calcium phosphates, a large number of works are devoted to the study of cationic substitutions in calcium phosphates, especially in hydroxyapatites, by DFT. DFT calculations make it possible to determine the most probable substitution position for impurity ions and monitor changes in the unit cell geometry, changes in the local positions of surrounding ions caused by impurities, etc. [20][21][22][23][24][25]. In most cases, DFT calculations help to confirm the assumptions that are made on the basis of experimental data. In addition to calculating positions and other geometric properties and cell geometry, DFT allows calculating various parameters that are associated with different physical properties of the substance. As an example, calculations of the electronic band structure, optical properties, energy of formation of defects, density of states, and EPR parameters, such as g-factor and hyperfine interaction constant A for HAp, are presented in [26][27][28][29]. Our work relates to the calculation of physical parameters that appear in the spectra of electron paramagnetic resonance.

Calculation of Electron Density
To calculate the electron density and configuration of paramagnetic impurities, the variational self-consistent field method was chosen. This method is implemented in the Quantum Espresso software package (PWscf module) using a plane waves (PWs) basis set and pseudopotentials to represent electron-ion interactions to the solution of the density functional theory (DFT) problem. The codes are constructed around the use of periodic boundary conditions, which allows for a straightforward treatment of infinite crystalline systems, and an efficient convergence to the thermodynamic limit for aperiodic but extended systems, such as liquids or amorphous materials [30].
The initial parameters of the crystal structure of hydroxyapatite were taken from the work [31]. Optimization was performed in two steps. In the first step, one hydroxyapatite cell was optimized. During the optimization of the geometry of one cell, the cell sizes and atomic positions were completely relaxed. In the second step, a supercell was obtained by doubling the new coordinates. Further, the calcium ion was replaced by a manganese ion in this structure. After this, the second optimization was performed with a fixed cell size. This was performed in order to take into account the locality of the substitution ion in further calculations, as the other surrounding cells of hydroxyapatite (in our calculations, there are about 20,000) do not change their size without impurity.
The unit cell was doubled along the A crystallographic axis (2 × 1 × 1). The size of the doubled cell is 18.848 × 9.424 × 6.879 Å. The impurity paramagnetic center replaces the Ca ion in one of two possible positions in the double cell according to the literature [32]. In our case, the Ca (I) position corresponds to the Ca (5) ion, and the Ca (II) position corresponds to the Ca (10) (Figure 1). Then, the cell geometry is optimized by energy minimization. This approach allows obtaining configurations of impurity centers in a crystal cell and numerical values of the electron density distribution. The positions of the ions could be varied, and the volume and size of the original cell were preserved.
The unit cell was doubled along the A crystallographic axis (2 × 1 × 1). The size of the doubled cell is 18.848 × 9.424 × 6.879 Å. The impurity paramagnetic center replaces the Ca ion in one of two possible positions in the double cell according to the literature [32]. In our case, the Ca (I) position corresponds to the Ca (5) ion, and the Ca (II) position corresponds to the Ca (10) (Figure 1). Then, the cell geometry is optimized by energy minimization. This approach allows obtaining configurations of impurity centers in a crystal cell and numerical values of the electron density distribution. The positions of the ions could be varied, and the volume and size of the original cell were preserved.
Calculations were implemented using plane-wave pseudopotential using the Purdue-Burke-Ernzerhof (PBE) exchange-correlation functional [33]. Structural optimization was performed with the kinetic energy cutoffs of 60 Ry for the smooth part of the electron wave functions and 400 Ry for the augmented electron density (as defined by the convergence test) and with a 1 × 1 × 1 Monkhorst-Pack k-point mesh.
The electron density distribution is shown in Figure 1. The figure shows that the electron density in the region of the Mn 2+ ion is spherical, but the electron density of the nearest PO4 groups is far from spherical. Therefore, the use of the model of point charges, in this case, is unjustified.

Calculation of the Parameters of the Electric Crystal Field
In the point-charge model, the crystal field includes the electric fields from all charged ions in the environment. The crystal field parameters can be calculated through the Coulomb interaction: where = [34]. As described above, it is insufficient to use the model of point charges for the HAp+Mn 2+ system. Therefore, this model was modified. The interaction of the impurity ion with the cores of the ions and with the electron density was taken into account separately: where j runs over all the ion cores in the lattice, ZCa = 10, ZO = 6, ZP = 5, and ZH = 1. The integral is taken over the entire electron density. Calculations were implemented using plane-wave pseudopotential using the Purdue-Burke-Ernzerhof (PBE) exchange-correlation functional [33]. Structural optimization was performed with the kinetic energy cutoffs of 60 Ry for the smooth part of the electron wave functions and 400 Ry for the augmented electron density (as defined by the convergence test) and with a 1 × 1 × 1 Monkhorst-Pack k-point mesh.
The electron density distribution is shown in Figure 1. The figure shows that the electron density in the region of the Mn 2+ ion is spherical, but the electron density of the nearest PO 4 groups is far from spherical. Therefore, the use of the model of point charges, in this case, is unjustified.

Calculation of the Parameters of the Electric Crystal Field
In the point-charge model, the crystal field includes the electric fields from all charged ions in the environment. The crystal field parameters B q k can be calculated through the Coulomb interaction: where r k i = ∞ 0 r k R 2 3d r 2 dr [34]. As described above, it is insufficient to use the model of point charges for the HAp+Mn 2+ system. Therefore, this model was modified. The interaction of the impurity ion with the cores of the ions and with the electron density was taken into account separately: where j runs over all the ion cores in the lattice, Z Ca = 10, Z O = 6, Z P = 5, and Z H = 1. The integral is taken over the entire electron density. The calculation of the crystal field parameters was carried out using a computer program written in the Python language for a sphere of double cells, in the center of which there is an impurity ion. The calculated positions of the impurity center for two possible variants of substitution at the positions Ca (I) and Ca (II) were taken as the origin of coordinates (centers of the spheres). This sphere includes entire double cells that are no further than a certain radius. Thus, proceeding from the electroneutrality of one cell, the electroneutrality of the entire crystal is taken into account. To maintain symmetry, the calculations use an odd number of cells in each direction. To accelerate the calculations, the following models are combined: the calculation of the parameters B k q taking into account the distributed (calculated) electron density inside the sphere of the nearest environment and the calculation of the parameters B k q in the point charges model for the remote environment. In the intermediate case, a light model of the electron density was used: the several points of the electron density are replaced by one point with average coordinates (reduction of sampling).
The calculation of electric crystal field parameters in the model with electron density was made in a radius from 0 to 6 nm (1727 double cells), an intermediate model was applied in a radius from 6 to 9 nm (3988 double cells), and the model of point charges was from 9 to 15 nm (21,064 double cells). The resulting parameters B q k were adjusted according to the article [35]. The results are shown in Table 1. To check the convergence, an increase in the radius of the spheres by 1.5 times was used; the change in the parameters was 0.01%. The calculation was also performed for different positions of Ca (I) and Ca (II) in the supercell; the change in the parameters did not exceed 1%.

HAp Samples and EPR Technique
Powders of hydroxyapatite doped with Mn 2+ with the chemical formula Ca 9 . 995 Mn 0 . 005 (PO 4 ) 6 (OH) 2 were synthesized by the wet precipitation technique at the Faculty of New Materials Science, Moscow State University. Details of the synthesis and powder characterization are given elsewhere [14,36]. The powders were investigated by EPR. The EPR spectra of impurity manganese ions in hydroxyapatite were obtained by exploiting the abilities of Bruker Elexsys 580/680 spectrometer in the X-(9 GHz) and W-(94 GHz) bands in the pulsed mode [9]. The EPR absorption spectra are recorded by the detection of the electron spin echo amplitude using the two-pulse Hahn sequence (π/2 ← 240 ns → π with the duration of the π/2 pulse equal to 24 ns) To acquire an electron spin echo signal with the appropriate signal-to-noise ratio, the EPR measurements were done at low temperatures by using a helium-flow cryostat.

Modeling of "Powder" EPR Spectra of Mn 2+
The Mn 2+ ion is in the 6 S 5/2 spin state and the fine structure in the EPR spectra should be absent in the first approximation. In [37], excited electronic states with nonzero spinorbital angular momentum are used to explain the appearance of ZFS. Using the results of this work in our case, we obtained a correction factor of the order 10 −7 , which is the same for all simulated spectra. Further, for the best agreement of the experimental data, this coefficient was chosen more accurately (for B 2 as 2 × 10 −7 , for B 4 as 0.7 × 10 −7 ). The obtained parameters of the crystal field were used to simulate the powder EPR spectra of manganese ions in hydroxyapatite in the X-and W-ranges. The modeling of the EPR spectra was carried out in the EasySpin module in MATLAB. The spectra are shown in Figure 2. The simulation was carried out simultaneously in the X and W bands of the EPR spectrometer, as this made it possible to isolate the ZFS effects. The simulation parameters of the experiment (field sweep, temperature, number of points) corresponded to the experimental ones. System parameters: g-factor from calculations in the Quantum Espresso program (for the Ca (I) center g x = 2.0082, g y = 2.0082, g z = 2.0078; for the Ca (II) center g x = 2.0096, g y = 2.008, g z = 2.0077), hyperfine interaction constant A from literature data (we consider it isotropic) A = 252 ± 3 MHz [14]. The intensity ratio of the manganese spectra at the Ca (I)/Ca (II) positions was 1.4. be absent in the first approximation. In [37], excited electronic states with nonzero spinorbital angular momentum are used to explain the appearance of ZFS. Using the results of this work in our case, we obtained a correction factor of the order 10 −7 , which is the same for all simulated spectra. Further, for the best agreement of the experimental data, this coefficient was chosen more accurately (for B2 as 2 × 10 −7 , for B4 as 0.7 × 10 −7 ). The obtained parameters of the crystal field were used to simulate the powder EPR spectra of manganese ions in hydroxyapatite in the X-and W-ranges. The modeling of the EPR spectra was carried out in the EasySpin module in MATLAB. The spectra are shown in Figure  2. The simulation was carried out simultaneously in the X and W bands of the EPR spectrometer, as this made it possible to isolate the ZFS effects. The simulation parameters of the experiment (field sweep, temperature, number of points) corresponded to the experimental ones. System parameters: g-factor from calculations in the Quantum Espresso program (for the Ca (I) center gx = 2.0082, gy = 2.0082, gz = 2.0078; for the Ca (II) center gx = 2.0096, gy = 2.008, gz = 2.0077), hyperfine interaction constant A from literature data (we consider it isotropic) A = 252 ± 3 MHz [14]. The intensity ratio of the manganese spectra at the Ca (I)/Ca (II) positions was 1.4.

Conclusions
It was demonstrated that the parameters of the crystal field in HAp can be calculated in the classical model using the distributed electron density. We tested the extraction for the two distinct calcium positions in HAp crystal field parameters to fit the EPR spectra of the doped Mn 2+ ions. Good agreement between the modeled and experimental data can be achieved by only suggesting that, in the hydroxyapatite powder under study, manganese ions replace calcium ions in both possible positions even at low concentration with the slight prevalence for the Ca (I) site. The obtained results show a way to extract information for the localization of paramagnetic ions in the HAp structure from the EPR spectra.
At present, there are different approaches to calculating the parameters of the crystal field. These approaches require complex mathematical calculations. The advantage of the presented in this paper technique is the simplicity of the mathematical model and the absence of complex mathematical calculations. To calculate the parameters of the crystal field, it is enough to set the initial positions of the ions in the system and indicate for which center the calculation will be carried out.

Conclusions
It was demonstrated that the parameters of the crystal field in HAp can be calculated in the classical model using the distributed electron density. We tested the extraction for the two distinct calcium positions in HAp crystal field parameters to fit the EPR spectra of the doped Mn 2+ ions. Good agreement between the modeled and experimental data can be achieved by only suggesting that, in the hydroxyapatite powder under study, manganese ions replace calcium ions in both possible positions even at low concentration with the slight prevalence for the Ca (I) site. The obtained results show a way to extract information for the localization of paramagnetic ions in the HAp structure from the EPR spectra.
At present, there are different approaches to calculating the parameters of the crystal field. These approaches require complex mathematical calculations. The advantage of the presented in this paper technique is the simplicity of the mathematical model and the absence of complex mathematical calculations. To calculate the parameters of the crystal field, it is enough to set the initial positions of the ions in the system and indicate for which center the calculation will be carried out.
The disadvantages of the method include the fact that the average radius is involved in the calculations. In fact, the radius of the impurity ion depends on the angle and time. In addition, the utilization of various software is required for the calculations, which implies high software costs and high qualifications of the researcher. Moreover, due to limitations in computer power, we cannot calculate the parameters of the crystal field at real (low) concentrations of dopants. On the other hand, our approach allows taking into account imperfections of real crystals, for example, by introducing the parameters of the distribution of the fine structure associated with dislocations in the crystal. The corresponding calculations are being prepared for publication.

Conflicts of Interest:
The authors declare no conflict of interest.