1 Introduction

Calcium phosphates have been widely studied for dentistry and orthopedic implants because of their high biocompatibility and bioresorbability [1–3]. Hydroxyapatite derivatives (Ca₅(PO₄)₃(OH)) (including carbonated apatitic phases) give bones, teeth and tendons their hardness and stability. By “playing” with the chemical composition, selective protein adsorption properties can be obtained for nanocrystalline hydroxyapatite [4]. Apatitic materials may also act as drug delivery systems [5–8].

Solid-state NMR is definitely suitable for the in-depth characterization of calcium phosphate phases. ³¹P NMR has been largely used in the past, in combination with ¹H → ³¹P Cross-Polarization transfer [9]. Moreover, inorganic components in bones were investigated [10]. ¹H solid-state NMR has been rarely used [11] due to the lack of resolution related to the strong homogeneous ¹H–¹H dipolar interaction. Multiple pulses experiments combined with MAS (CRAMPS) [12] have been implemented for the characterization of calcium phosphate phases [13]. Nevertheless, full assignment of proton lines remains difficult, though empirical trends correlating the isotropic chemical shifts and the H-bond networks have been proposed in the literature some decades ago [14]. In this paper, we demonstrate that first-principles calculations are validated for ³¹P solid-state NMR. Calculated CSA tensors are in very good agreement with experimental data obtained by MAS or static experiments. When considering ¹H MAS spectra (at very high field and very fast rotation frequency), calculated data become essential for full assignment of the experimental spectra. All calculations were done by using the GIPAW method first developed by Mauri and Pickard (see Section 2) [15].

¹⁷O solid-state NMR experiments appear as a straightforward tool of characterization of the local structure of inorganic materials [16] by probing directly and quantifying the different types of oxo-bridges present in the structure. The ¹⁷O isotropic chemical shift values of M–O–M′ fragments (M, M′ = Si, P, Ti, Zr,…) are indeed highly sensitive to the chemical nature of M and M′ [17–19]. However, the main drawback of this technique is its poor sensitivity, due to the low natural abundance of the ¹⁷O isotope (0.037%). In the field of phosphate compounds, this problem can nonetheless be overcome either by using ¹⁷O-enriched water for the hydrolysis step of the precursors [20] or by heating the final material under ¹⁷O-enriched water vapour [21] or ¹⁷O₂ gas [22]. Despite isotopic enrichment and high-resolution MQ-MAS (Multiple-Quantum Magic-Angle Spinning [23]) methods, the interpretation of ¹⁷O spectra still remains difficult, as no empirical rules have yet been found to enable the assignment of the peaks and to obtain structural information for each site. In this context, first-principles calculations can be particularly useful in the understanding of spectra, since accurate calculation of the NMR parameters will allow the assignment of all sites.

⁴³Ca (I = 7/2, 0.145%) NMR remains exotic in inorganic chemistry due to very low sensitivity and large quadrupole moment. It must be noticed that recent work dealing with model compounds and ⁴³Ca NMR characterization at high field has been carried out in the literature by using large amount of samples [24].

The outline of the article is the following: Section 2 presents the syntheses of the studied calcium phosphate phases, including protonated structures (CaHPO₄·2H₂O, Brushite and Ca(H₂PO₄)₂·H₂O, MCPM) and non-protonated structures (β-Ca(PO₃)₂ and Ca₄P₂O₉, TTCP). The first-principles calculations and the GIPAW approach are also described in Section 2. Section 3 deals with results and discussion: the calculated and experimental ³¹P are first discussed for all structures, including polymorphs of Ca(PO₃)₂. ¹H NMR data are discussed for protonated phases, emphasizing the calculation of CSA tensors and the strong influence of H bonding. Finally, ¹⁷O predicted data (isotropic chemical shifts and quadrupolar constants) are summarized and discussed.

2 Experimental

Brushite, CaHPO₄·2H₂O (Fig. 1), was synthesized by mixing calcium carbonate and H₃PO₄ (85%) in water (Ca/P = 1). The mixture was stirred at room temperature until precipitation. For Ca/P = 0.5, MCPM, Ca(H₂PO₄)₂·H₂O (Fig. 2), was obtained. By heating MCPM at 600 °C (10 h), β-Ca(PO₃)₂ (Fig. 3) was synthesized. TTCP, Ca₄P₂O₉ (Fig. 4), was obtained by mixing calcium carbonate and monetite (CaHPO₄) at 1400 °C for 10 h, followed by quenching in air. The obtained compounds were systematically checked by powder X-ray diffraction (Brushite, JCPDS: 72-1240; MCPM, JCPDS: 70-0090; β-Ca(PO₃)₂, JCPDS: 11-0039; TTCP, JCPDS: 70-1379).

¹H NMR experiments were recorded on a Bruker AVANCE 750 spectrometer in order to enhance resolution (¹H: 750.16 MHz). A zirconia rotor (2.5 mm) was spun at 33 kHz at the magic angle. The number of scans (NS) was 4, the pulse length was 1.7 μs (30°), and the relaxation delay (RD) was 30 s. The spectra were referenced to TMS. ³¹P MAS NMR experiments were performed on a Bruker AVANCE 300 (³¹P: 121.44 MHz) using a 4-mm Bruker probe under ¹H high power decoupling (TPPM pulse sequence) (pulse length: 1.6 μs (30°), RD: 30 s). Spectra were referenced to H₃PO₄ (85%). A zirconia rotor (4 mm) was spun at 3 and 14 kHz.

The first-principles calculations were performed within the Kohn–Sham DFT. The PBE generalized gradient approximation [25] was used and the valence electrons were described by norm-conserving pseudopotentials [26] in the Kleinman–Bylander [27] form. The core definition for O is 1s² and it is 1s²2s²2p⁶ for P and Ca. The core radii are 1.2 a.u. for H, 1.5 a.u. for O and 2.0 a.u. for P and Ca. The wave functions are expanded on a plane wave basis set with a kinetic energy cut-off of 80 Ry. The crystalline structure is described as an infinite periodic system using periodic boundary conditions. The NMR calculations were performed for the experimental geometries determined by neutron or X-ray diffraction for the various calcium phosphates: brushite, CaHPO₄·2H₂O [28]; MCPM, Ca(H₂PO₄)₂·H₂O [29]; β-Ca(PO₃)₂ [30]; γ-Ca(PO₃)₂ [31]; Ca₄P₂O₉ [32]. In the case of Ca₄P₂O₉, a significant better agreement is obtained between experimental and calculated ³¹P parameters by taking into account a geometry optimization (starting from the experimental geometry and allowing the positions of the atoms to relax using DFT calculations). The new coordinates obtained after relaxation are summarized in Supporting information. The integral over the first Brillouin zone was performed using a Monkhorst–Pack 2 × 2 × 2 k-point grid [33] for the charge density and electric-field gradient calculation and a 4 × 4 × 4 k-point grid for the chemical shift tensor calculation. The calculations have been performed at the IDRIS supercomputer centre of the CNRS using a parallel IBM Power4 (1.3 GHz) computer; the calculation of the Electric Field Gradient (EFG) and CSA tensors requires from 2 to 10 h on 32 processors. The shielding tensor is computed using the GIPAW approach which permits the reproduction of the results of a fully converged all-electron calculation, while the EFG tensors are computed using a PAW approach. The isotropic chemical shift δ_iso is defined as δ_iso = −[σ − σ^ref] where σ is the isotropic shielding and σ^ref is the isotropic shielding of the same nucleus in a reference system. In our calculations, absolute shielding tensors are obtained. The σ^ref values were chosen by comparison between experimental and calculated δ_iso in a series of SiO₂ polymorphs [34a] for ¹⁷O and other organophosphorus derivatives for ¹H and ³¹P: σ^ref (¹⁷O) = 261.5 ppm, σ^ref (¹H) = 31.0 ppm, and σ^ref (³¹P) = 300.7 ppm. Diagonalisation of the symmetrical part of the calculated shielding tensor provides its principle components δ₁₁, δ₂₂ and δ₃₃ defined such as |δ₃₃ − δ_iso| ≥ |δ₁₁ − δ_iso| ≥ |δ₂₂ − δ_iso| and δ_iso = 1/3(δ₁₁ + δ₂₂ + δ₃₃). Similarly, the principle components V_xx, V_yy, and V_zz of the EFG tensor defined as |V_zz| ≥ |V_xx| ≥ |V_yy| are obtained by diagonalisation. The quadrupolar interaction can then be characterized by the quadrupolar coupling constant C_Q (with its sign) and the asymmetry parameter η_Q defined as: C_Q = eQV_zz/h and η_Q = (V_yy − V_xx)/V_zz. It has to be noticed that the role of Ca–O bonds within the frame of GIPAW calculations has been carefully studied by Profeta et al. [34b]. It was demonstrated that an easy and transferable correction of the PBE–DFT theory was possible for the accurate calculations of NMR parameters.

3 Results and discussion

3.2 ¹H: calculated vs. experimental data

Ab initio calculations of ¹H parameters are rarely described [38]. In most cases, ¹H data are calculated within the cluster approximation [39], with some exceptions involving mainly organic molecules [40]. Only δ_iso (¹H) values are discussed. It should be noticed that empirical correlations involving the δ_iso (¹H) values and the internuclear distances have been previously described in the literature [14]. The goal of this section is to revisit such correlations from the first-principles calculation point of view.

Fig. 1b presents the ¹H MAS NMR spectrum of brushite obtained at very high field (B₀ = 17.6 T) and very fast rotation frequency (up to 34 kHz). Full resolution is not attained due to the strong ¹H/¹H dipolar interaction, which is homogeneous in nature [41]. However, three resonances located at ∼4.0, 6.5 and 10.0 ppm are observed. Calculated GIPAW values are summarized in Table 1, leading to the calculated spectrum shown in Fig. 1b. The most deshielded line corresponds obviously to POH groups. The other resonances are safely assigned to protons belonging to water molecules. In a sense, calculations allowed full assignment and interpretation of the MAS spectrum. Calculated ¹H data related to MCPM are also given in Table 1.

Fig. 5 presents all calculated δ_iso (¹H) values vs. the shortest OH⋯O distances (extracted mainly from neutron diffraction data). Distinct ranges for H₂O and POH resonances are evidenced. Moreover, for a given range, deshielding of protons associated with a shortening of the OH⋯O distances is clearly demonstrated. Recently [42], such results have been extended to ¹H data related to phenylphosphonic and phenylphosphinic acids. The ¹H CSA data shown in Table 1 are consistent with nearly axial symmetry for all tensors (H₂O and POH). Such data are in agreement with results already published in the literature [43]. Nevertheless, Wu et al. [44] have pointed out that ¹H CSA tensors could deviate from axial symmetry in the case of crystalline hydrates (involving perchlorate and sulphate derivatives). Fig. 6 presents the absolute orientations of ¹H CSA tensors in the case of MCPM. For POH and H₂O protons, the unique axis (δ₃₃) of the CSA tensors is nearly collinear to the OH⋯O bond direction.

3.3 ¹⁷O: calculated data

In the field of calcium phosphates, few ¹⁷O NMR results have been published until now to our knowledge. Only ¹⁷O MAS NMR spectra of Brushite and hydroxyapatite were previously reported [45].

The aim of the calculation of ¹⁷O NMR parameters in the various compounds presented in this paper is to determine whether there is a possibility to distinguish clearly the various types of oxygen environments present in calcium phosphates (Ca⋯O, P–O⋯Ca, P–O–P, P–OH and H₂O) by using ¹⁷O NMR spectroscopy. The results are summarized in Table 1. The previously reported ¹⁷O NMR values for P–O–P sites in h-P₂O₅ [22], Ca⋯O in lime [46] (CaO) and P–O⋯Ca environments in brushite and HAp45 are in very good agreement with the ones obtained by the calculations. Moreover, it should be noticed that the relatively large ¹⁷O quadrupolar coupling constants calculated for water molecules in brushite and MCPM compare well with those experimentally observed in p-toluenesulphonic acid hydrate [47] and oxalic acid dihydrate [48]. Similarly, the ¹⁷O P–OH parameters are consistent with those measured in phenylphosphonic [49] and phenylphosphinic [42] acid.

Ca⋯O environments are expected in a very distinct area with high chemical shift values (around 300 ppm) and very small quadrupolar coupling constants (less than 1 MHz). Similarly H₂O sites are well separated with small chemical shift values (≤0 ppm) and very large quadrupolar coupling constants (>7.5 MHz). On the other hand, P–O⋯Ca, P–O–P and P–OH environments appear in the same chemical shift range (between 65 and 115 ppm). To evaluate which type of oxygen environment can be clearly distinguished experimentally, we have calculated the ¹⁷O 3QMAS spectrum corresponding to the different oxygen sites (Fig. 7) [50]. Each ¹⁷O site gives rise to a sharp resonance in the δ_3Qiso dimension at the position:

δ3Qiso=−17/31δiso+10/31δiso2Q

and a broad MAS spectrum in the MAS dimension with a centre of gravity given by:

δ3QMAS=δiso+δiso2Q

where δ_iso is the isotropic chemical shift and δiso2Q the second-order quadrupolar shift given for I = 5/2 by 6000PQ2ω02, with ω₀ the Larmor frequency and PQ=CQ(1+ηQ2/3)1/2.

The 2D spectrum shows that all types of oxygen environments investigated in this study appear in relatively distinct areas and should therefore be distinguishable using the 3QMAS NMR sequence. Indeed, even if P–O⋯Ca, P–O–P and P–OH environments have similar isotropic ¹⁷O chemical shifts, the discrimination is possible, thanks to different C_Q values.

It can be noticed that the ¹⁷O chemical shift range observed for P–O⋯Ca environments is rather large and it was suggested that the shifts could discriminate between different P–O bond lengths in PO₄³⁻ units [45]. Calculations show indeed that δ_iso (¹⁷O) tends to increase with P–O⋯Ca bond length.

4 Conclusions

It has been shown that first-principles calculations (using the GIPAW approach) are suitable for the interpretation of multinuclear solid-state NMR spectra of calcium phosphate phases, even in the case of non-fully resolved ¹H MAS experiments. Full tensor data are obtained for ¹H and ³¹P, including the absolute orientation of the principle axes. Interesting correlations between the δ_iso (¹H) values and the shortest OH⋯O distances have been established, giving a more theoretical background to empirical data already published in the literature. Calculated ¹⁷O parameters tend to prove that the various oxygen sites can be separated in MQ-MAS experiments. Work is in progress in the laboratory for recording ¹⁷O spectra of enriched calcium phosphate phases.

Acknowledgments

Dr. D. Massiot (CRMHT) is warmly acknowledged for helpful discussions and for the access to the 750-MHz spectrometer in Orléans, France. The calculations have been performed in the IDRIS supercomputer centre of the CNRS (Project 51461).

Appendix Supporting information

The table includes fractional atomic coordinates of Ca₄P₂O₉ structure after relaxation. The lattice parameters were not relaxed; a = 7.023 Å, b = 11.986 Å, c = 9.473 Å. The symmetry of the crystal was also constrained to be monoclinic P1211.