Molecular Dynamics Simulations of Hydroxyapatite Nanopores in Contact with Electrolyte Solutions: The Effect of Nanoconfinement and Solvated Ions on the Surface Reactivity and the Structural, Dynamical, and Vibrational Properties of Water

Abstract

Hydroxyapatite, the main mineral phase of mammalian tooth enamel and bone, grows within nanoconfined environments and in contact with aqueous solutions that are rich in ions. Hydroxyapatite nanopores of different pore sizes (20 Å ≤ H ≤ 110 Å, where H is the size of the nanopore) in contact with liquid water and aqueous electrolyte solutions (CaCl₂ (aq) and CaF₂ (aq)) were investigated using molecular dynamics simulations to quantify the effect of nanoconfinement and solvated ions on the surface reactivity and the structural and dynamical properties of water. The combined effect of solution composition and nanoconfinement significantly slows the self-diffusion coefficient of water molecules compared with bulk liquid. Analysis of the pair and angular distribution functions, distribution of hydrogen bonds, velocity autocorrelation functions, and power spectra of water shows that solution composition and nanoconfinement in particular enhance the rigidity of the water hydrogen bonding network. Calculation of the water exchange events in the coordination of calcium ions reveals that the dynamics of water molecules at the HAP–solution interface decreases substantially with the degree of confinement. Ions in solution also reduce the water dynamics at the surface calcium sites. Together, these changes in the properties of water impart an overall rigidifying effect on the solvent network and reduce the reactivity at the hydroxyapatite-solution interface. Since the process of surface-cation-dehydration governs the kinetics of the reactions occurring at mineral surfaces, such as adsorption and crystal growth, this work shows how nanoconfinement and solvation environment influence the molecular-level events surrounding the crystallization of hydroxyapatite.

1. Introduction

Crystallization is a phenomenon that is strongly dependent on the environment in which it occurs. Specifically, it is well-recognized that many features of crystal growth are affected by confinement [1,2,3,4,5,6,7] and solution composition [2,8,9,10,11,12]. Confinement can enable control over size, polymorphism and orientation, single crystal/polycrystalline structure, and morphology [13]. At the same time, modification of the solvation environment through the addition to the solution of organics [14], peptides [15], polymers [16], or simple ionic salts [17,18,19,20], can influence the kinetics and thermodynamics of crystal nucleation and growth. A perfect example of a crystallization process that occurs within a confined environment and in multicomponent aqueous solutions is biomineralization [2,21]. In particular, hydroxyapatite (HAP), molecular unit formula (Ca₁₀(PO₄)₆(OH)₂), is arguably the most important biomineral as it represents the main constituent of bone [22].

It is commonly accepted that bone can be described by three main levels of porosity [23], which are nested hierarchically one inside another as a set of Russian dolls in microcirculatory pathways [24]. The macroscopic porous network corresponds to the vascular (or Haversian) porosity, which consists of the Havers and Volkmann canals (typical diameter of 50 µm). The mesoscale is the lacunar-canalicular porosity made of the osteocytic lacunae and canaliculi channels. This porosity contains the osteocytes’ stellar network of bone (typical pore size of 100 nm) and the smallest porosity level in bone corresponds to the spaces inside the collagen-apatite structure (typical pore size of 5 nm). At this nanoconfined environment, water strongly interacts with the ionic crystal [25,26] and plays a key role in structuring the apatite mineral [27,28] and during mineralization [29] by acting as a prominent charge carrier, maintaining the pH of the medium, and transporting ions [26].

Determining the structure, dynamics, and reactivity of fluids in HAP nanopores is therefore of great interest to improve our understanding of bone. Molecular dynamics (MD) simulations of mineral-water interfaces can provide a unique probe of the processes controlling crystallization [27], such as the water exchange around the constituent cations (Ca²⁺ in hydroxyapatite), which is generally accepted to be the rate limiting step to crystal growth from aqueous solution of ionic crystals [28,29,30].

Previous computational investigations of HAP focused on its bulk and surface properties [31,32,33] as well as on the early stages of homogeneous nucleation from solution [34,35], but the collective effect of nanoconfinement and solution chemistry on the molecular-level processes controlling the surface reactivity of HAP has been entirely overlooked by molecular modelling. This aspect is particularly relevant to further our understanding of bone formation, which begins with the growth of ultrathin HAP platelets within confined volumes in the collagen fibrils [36,37] and in aqueous environments that are far from pure water but rich in ions.

This study reports MD simulations of HAP nanopores of different pore sizes (20 Å ≤ H ≤ 110 Å) in contact with liquid water and aqueous electrolyte solutions (CaCl₂ (aq) and CaF₂ (aq)). The aim is to quantify the effect of confinement and solution chemistry on the frequency of hydration–dehydration around the calcium surface sites as well as on the structural, dynamical, and vibrational properties of water. The realism of this simulation set-up represents a significant improvement over previous computational work on hydrated HAP.

2. Computational Methods

2.1. Simulation Details

MD simulations were performed using version 4.03.3 of the DL_POLY computational chemistry package [38]. The velocity Verlet algorithm with a time step of 0.1 fs was used to integrate the equations of motion. The Melchionna modification of the Nosé-Hoover algorithm [39] was used to maintain a temperature of 300 K and a pressure of 105 Pa (1 bar) with 0.5 ps for the thermostat (NPT ensemble) and barostat relaxation times. The electrostatic forces were calculated by means of the Smooth Particle Mesh Ewald method [40]. A 9.0 Å cut-off was used for the real-space part of the Ewald sum and the Van der Waals forces. Periodic boundary conditions were applied in all three directions of the unit cell.

2.2. Interatomic Potential Model

A combination of potential models was employed to describe the interatomic interactions in this multi-component system. The interatomic potential model developed by de Leeuw for modelling apatite crystals [41] was used for HAP. In this forcefield, phosphate and hydroxyl group bonds are described as the sum of a Morse and a Coulombic potential, phosphate bond angles by a harmonic potential, and non-bonded interactions by Buckingham potentials. This forcefield makes use of a shell model in which each oxygen anion in the phosphate and hydroxyl groups consists of both a core and a massless shell connected by a spring, in order to model the atom’s electronic polarizability. The water molecules were described using the simple point charge (SPC/E) model [42] because it gives a density, radial distribution functions, self-diffusion coefficient for water, and hydrogen-bond dynamics in good agreement with the experiment [42,43]. HAP–water interactions were described using non-bonded potentials derived in previous MD studies [44,45], which are compatible with the current HAP shell model. The interactions between Ca²⁺–X⁻ and X⁻–X⁻ pairs in solution (where X⁻ designates Cl⁻ or F⁻), and between the ions in solution and the phosphate and hydroxyl groups of HAP were described using the Buckingham potentials developed by Rabone and de Leeuw for modeling natural apatite crystals [46]. Finally, the halide–water interactions were described using the force field parameters derived by Dang [47], because it provides an accurate description of the diffusion coefficients of halide ions, and of the dynamical and structural properties of the hydrogen bonding network [48,49,50,51,52]. The list of parameters and an assessment of the forcefield used in this work are reported in the Supplementary Materials, Sections S1 and S2.

2.3. Simulation Protocol

Hydroxyapatite, (Ca₁₀(PO₄)₆(OH)₂), is viewed as a hexagonal primitive cell with P6₃/m space group. The initial coordinates and crystal cell parameters were taken from Sudarsanan and Young [53] and a 3 × 3 × 4 supercell was then generated to conduct the simulations. HAP nanopores with initial pore sizes, denoted by H, equal to 20, 40, 60, and 110 Å were generated by varying the c-axis of the crystal, with the resulting surface which corresponds to the (0001) basal plane. This model corresponds, to a good approximation, to a pore because at the atomic scale the pore surface will be relatively flat. More complex surfaces will be considered in future work. The vacuum was filled with water corresponding to the experimental density of 1 g·mol⁻¹. All generated systems were equilibrated for 50 ps in the microcanonical (NVE) ensemble, followed by 100 ps in the NPT ensemble, where the volume of the HAP nanopores was monitored to verify its convergence and confirm that the system was at equilibrium. Finally, the production phase in the NPT ensemble lasted for 3 to 4 ns. Starting from the last configuration of the HAP (H = 60 Å)–water system, aqueous CaX₂ (X = F and Cl) solutions with concentrations ranging from 0 to 1 mol·kg⁻¹ were generated by randomly replacing N water molecules with N/3 calcium ions and 2N/3 anions. For each HAP-solution system, 0.5 ns of equilibration were followed by 3 to 4 ns of production (NPT ensemble). Figure 1 displays the equilibrated structure of the HAP (H = 60 Å)–water system. SPC/E bulk water and aqueous CaCl₂ solutions were also considered. The number of atoms and concentrations of the systems considered in the present study are listed in

Table 1. Bulk solutions and hydroxyapatatite nanopores considered in this study. Concentration of calcium chloride and calcium fluoride in mol·kg⁻¹ (b).

	Number of Ca10(PO4)6(OH)2 Units	Number of CaX2	Number of H2O	b
Solution
Water	–	–	729	0.0
CaCl2	–	8	705	0.6
	–	16	681	1.3
CaF2	–	8	705	0.6
	–	16	681	1.3
HAP, H = 110 Å
Water	36	–	2505	0.0
HAP, H = 60 Å
Water	36	–	1363	0.0
CaCl2	36	4	1351	0.2
	36	11	1330	0.5
	36	23	1294	1.0
CaF2	36	4	1351	0.2
	36	11	1330	0.5
	36	23	1294	1.0
HAP, H = 40 Å
Water	36	–	910	0.0
HAP, H = 20 Å
Water	36	–	455	0.0

.

2.4. Analysis of the Simulation Data

Interatomic bonding pairs were analyzed through the generation of radial distribution functions (RDF), g_αβ(r), which represent the probability, relative to a random distribution, of finding atoms of types α and β separated by a distance r [54].

The number of hydrogen bonds (HBs) between water molecules in bulk solutions and within HAP nanopores was extracted from the MD trajectories. In particular, in order to determine the existence of an HB between two water molecules, the following three geometrical criteria were applied [55]: (i) the oxygen distance is less than 3.5 Å; (ii) the intermolecular hydrogen–oxygen distance is less than 2.45 Å; (iii) the oxygen–oxygen–hydrogen angle is less than 30°.

The self-diffusion scalar coefficient of water molecules (D_W) was calculated using the mean square displacement (MSD) utility of the DL_POLY code, which accurately determines diffusion coefficients by processing the trajectory file generated and computing the mean square displacement for the different atomic species in the simulation using multiple time origins [38,54].

Vibrational properties of the water molecules relied on the calculation of the velocity-autocorrelation function (VACF) of oxygen and hydrogen atoms. VACF is defined as follows:

$$VACF\left(t\right)=\frac{1}{N_O{N}_{atm}}{\displaystyle \sum }_{j=1}^{N_O}{\displaystyle \sum }_{i=1}^{N_{atm}}{\mathit{v}}_i\left(t_j\right)\times {\mathit{v}}_i\left(t_j+t\right)$$

(1)

where $v_i$ is the velocity vector of O or H atoms in the water molecule $i$, $N_O$ is the number of time origins spaced by t, and $N_{atm}$ is the number of O or H atoms. From the Fourier transform of the VACF of each atom, the vibrational density of states (VDOS) of each atom was then calculated [54].

The dynamics of the first hydration shell of calcium ions in solution and in the HAP surface were characterized using the “direct” method proposed by Hofer and co-workers [56]. This method was successfully applied for, among other purposes, the characterization of the dynamics of the coordination shell of hydrated alkaline earth metal ions and their carbonate and bicarbonate complexes [57], as well as the quantification of the water exchange frequency at structurally distinct calcium sites on the calcite surface [58]. The MD trajectories were analyzed for water molecule movements and whenever a water molecule crossed the boundary of the calcium coordination shell, its path was followed; if its new position outside or inside this shell lasted for more than 0.5 ps, the event was counted as a real exchange event [56]. For calcium ions, the first shell was defined to fall within the first minimum of the calcium versus water oxygen (Ca–O_W) RDF (3.3 Å).

The speciation of Ca²⁺ and X⁻ in solution was determined using the following ion-pairing criteria: (i) a contact ion pair (CIP) when Ca²⁺ and X⁻ are in direct physical contact (r_Ca-Cl < 3.5 Å and r_Ca-F < 3.0 Å); (ii) a solvent-shared ion pair (SSHIP) when Ca²⁺ and X⁻ (X = Cl or F) are separated by one water molecule; (iii) a solvent separated ion pair (SSIP) when Ca²⁺ and X⁻ are separated by two water molecules; and (iv) free ions (FIP) when the above conditions are not met and the cation and anion are fully hydrated beyond the second hydration shell.

3. Results and Discussion

3.1. Structural, Dynamical, and Vibrational Properties of Water

Figure 2 reports the effect of confinement and salt concentration on the structure of water, which was quantified in terms of the water oxygen–oxygen (O–O) radial distribution function (RDF). In Figure 2a, the location of the first peak, which corresponds to the average position of two H-bonded water molecules, is not significantly affected by the level of water confinement in the HAP nanopores. However, as the size of the nanopore decreases, the intensity of the first peak becomes more pronounced and the position of the second maximum moves to lower distances. This behavior indicates structuring of water in the nanopores compared with bulk liquid. On the other hand, the presence of electrolytes in the solutions in contact with HAP nanopores (Figure 2b,c) and in bulk solutions (Figure 2d,e) decrease the first minima and increase the first and second maxima in the O–O RDFs. CaCl₂ and CaF₂ salts therefore disrupt the structure of water. In particular, the effect on the O–O structure is more pronounced in CaF₂ (aq), where for the 1.0 mol·kg⁻¹ solution the second shell has collapsed into the first one (Figure 2c). A similar behavior in the interatomic distances of water was observed from MD simulations of NaCl and NaF aqueous solutions [45]. The stronger electrostatic field of F⁻ compared with Cl⁻ can be used to explain the more pronounced effect of fluorine ions on the structure of water.

Further insights into the combined effect of confinement and solution composition can be determined in Figure 3 from the angular distribution functions (ADF), denoted by ρ(θ), of the O–O–O trios between three water molecules connected by hydrogen bonding. Each water molecule can form two HBs involving their own hydrogen atoms and two further HBs by utilizing the hydrogen atoms attached to neighboring water molecules (Figure 3). These four HBs optimally arrange themselves tetrahedrally (O–O–O angle of 109.5°) around each water molecule as found in ordinary ice. In pure water, the function ρ(θ) is centered at 109.5°, which is an indication of the local tetrahedral clustering. Another indication of water structuring in HAP nanopores is the narrowing of the angular distribution around 109.5°, especially for H = 20 Å and 40 Å. The presence in solution of calcium and chloride ions disrupts the orientational order in the hydrogen-bonding network and shifts the distribution of the O–O–O angle to lower values.

In order to rationalize the impact of nanopore confinement and salt concentration on the arrangement of the water molecules, the MD trajectories were analyzed to determine the distribution of the number of HBs.

Table 2. The distribution of the number of hydrogen bonds per water molecule. Results obtained from MD simulations of bulk water, aqueous CaCl₂ solutions, water within HAP nanopores of different sizes (20 Å ≤ H ≤ 110 Å), and CaCl₂ and CaF₂ solutions within the HAP nanopores with H = 60 Å. The values given are percentages of molecules with the given number of hydrogen bonds. Concentration in mol·kg⁻¹ (b).

		Number of Hydrogen Bonds (%)
	b	0	1	2	3	4	5	Average
Bulk solution
Water	–	0.0	0.9	8.8	33.0	51.3	5.9	3.53
CaCl2	0.6	0.2	2.7	14.5	34.0	43.7	4.8	3.30
	1.3	0.4	4.7	19.2	34.0	37.6	4.0	3.16
CaF2	0.6	0.2	2.8	13.5	33.4	44.7	5.3	3.36
	1.3	0.5	4.7	17.1	33.1	39.6	4.8	3.21
HAP, H = 110 Å
Water	–	0.4	2.8	11.1	32.8	47.4	5.5	3.41
HAP, H = 60 Å
Water	–	0.5	3.9	8.4	20.6	61.9	4.7	3.54
CaCl2	0.2	0.7	4.5	10.0	22.0	58.4	4.5	3.46
	0.5	0.6	5.0	13.3	24.8	52.1	4.0	3.35
	1.0	0.9	6.8	16.8	26.9	45.2	3.4	3.19
CaF2	0.2	0.9	5.3	11.1	27.1	50.4	5.1	3.36
	0.5	1.1	6.1	15.2	30.6	41.9	5.0	3.22
	1.0	1.7	8.6	17.2	27.3	40.8	4.2	3.10
HAP, H = 40 Å
Water	–	0.9	6.7	8.5	19.4	60.0	4.3	3.44
HAP, H = 20 Å
Water	–	0.8	11.2	12.4	12.4	59.7	3.5	3.29

shows that in bulk SPC/E water, the average number of HBs per water molecule is 3.53. This result is in very good agreement with the experimental value of 3.58 obtained from neutron diffraction data [59]. Confinement of water and, in particular, the presence of electrolytes reduce the average HBs per water molecule. The percentage of molecules that form four HBs decreases from 51% in pure water to 44% in the 0.6 mol·kg⁻¹ CaCl₂ solution, and to 38% in the 1.3 mol·kg⁻¹ CaCl₂ solution. Moreover, the number of molecules that form only two HBs is double of that in pure water (see Table 2). Similar behavior is observed for the CaF₂ solutions.

Nanoconfinement also affects the water H-bonding network. In the aqueous solutions within the HAP nanopores, water molecules forming three HBs decrease, and those involved in only two HBs substantially increase compared with bulk water (see Table 2), which explains the deviation from ideal tetrahedrality (109.5°) in the O–O–O ADF profiles reported in Figure 3. In particular, the average number of HBs decreases to only 3.3 per water molecule in the HAP (H = 20 Å) nanopores. Notice also the average value of 3.10 in the 1.9 mol·kg⁻¹ CaF₂ solutions within the HAP (H = 60 Å) nanopore, which is due to a substantial decrease of water molecules forming four HBs (40.8% compared with 51.3% in bulk water) and an increase of molecules forming only two (17.2%) and three (27.3%) HBs. This explains the broadening in CaF₂ solutions of the O–O–O ADF in CaF₂ (Figure 3c).

The self-diffusion coeffiecients of oxygen water, denoted by D_W, were used to quantify the effect of salt composition and nanoconfinement on the dynamics of water (

Table 3. The self-diffusion coefficient (D_W) for bulk and confined water. Concentration of CaCl₂ and CaF₂ in mol·kg⁻¹ (b).

	b	D W × 10−10 (m2·s−1)
Solution
Water	–	25.8
CaCl2	0.6	21.1
	1.3	16.3
CaF2	0.6	22.3
	1.3	17.1
HAP, H = 110 Å
Water	–	21.9
HAP, H = 60 Å
Water	–	5.4
CaCl2	0.2	3.9
	0.5	3.5
	1.0	2.9
CaF2	0.2	6.6
	0.5	5.3
	1.0	3.9
HAP, H = 40 Å
Water	–	2.9
HAP, H = 20 Å
Water	–	0.2

). The value of D_W for bulk SPC/E water is 25.8 × 10⁻¹⁰ m²·s⁻¹, in good agreement with previous work [60]. In bulk CaCl₂ and CaF₂ solutions there is a marked decrease in water diffusivity, which is substantially higher than in aqueous alkali halide solutions [61]. For example, in 1.9 mol·kg⁻¹ NaCl (aq) the value of D_W is 22.0 × 10⁻¹⁰ m²·s⁻¹, which is higher than the value of D_W (16.3 × 10⁻¹⁰ m²·s⁻¹) in the 1.3 mol·kg⁻¹ CaCl₂ (aq). This result can be rationalized in terms of the stronger coordination of water molecules to Ca²⁺ compared with Na⁺. In fact, the mean residence time (MRT) of water molecules in the first hydration shell of Ca²⁺ is in the range of 23–105 ps (depending on the interatomic potential model used to describe the calcium-water interaction [57,62,63]), whereas the MRT of Na⁺ is only 8 ps [61]. This implies that the hydration shell of the sodium ion is more labile than that of calcium ions. Previous computational studies showed that in nanoconfined environments the diffusion of water is considerably slower than in the bulk [63,64]. However, compared with simulations of graphene oxide layers [63], where the water diffusion in a layer spacing of 15 Å converged to that of bulk water, in the HAP (H = 110 Å)–water system the value of D_W (19.3 × 10⁻¹⁰ m²·s⁻¹) is still 15% lower than in the bulk water (see Table 3). In the HAP (H = 40 Å) and HAP (60 Å) systems, which have sizes more relevant to the nanopores found in the osteocytes’ stellar network of bone and in the spaces inside the collagen-apatite [25], the self-diffusion coefficients are 5.4 × 10⁻¹⁰ and 2.9 × 10⁻¹⁰ m²·s⁻¹, respectively, which are about 80% lower than in bulk water. The combined effect of electrolytes and nanoconfinement in the HAP (H = 60 Å)–aqueous CaX₂ solutions (X = Cl and F) systems further decreases the dynamics of water (see Table 3).

Further insights into the dynamics of water are obtained from the analysis of the velocity-autocorrelation function (VACF) of the water oxygen atoms. The occurrence of a minimum in the VACF profile results from the so-called “cage effect” for the tagged particle—i.e., it takes some time for the particle to escape from the cage formed by its surrounding neighbors [65]. The oscillatory behavior and position of the first minimum in the VACF profile can therefore be used to probe the interaction of the tagged particle with the surrounding cage [66,67]. In the VACF profile of bulk water (Figure 4a), the first minimum located at 85 fs is followed by an oscillation attributed to the intermolecular O–O stretch vibration of two water molecules linked by H-bonding [68]. Changes in the profiles of the VACF of water with solution composition and nanoconfinement can be used as descriptors of water mobility. Nanoconfinement in particular leads to a marked oscillatory behavior of the VACF (Figure 4b), which is indicative of strengthening of the HB water network on going from H = 110 Å to H = 20 Å. The presence in solution of calcium and chloride ions causes a similar but less pronounced effect on the VACF profiles (Figure 4c).

Figure 5 reports the simulated total vibrational density of states (VDOS) of water, which was computed as the sum of the Fourier transform of the VACF of the oxygen and hydrogen atoms. At low frequencies, the spectrum of liquid water has two characteristic peaks resulting from the H-bonding network, one centered at 50 cm⁻¹, corresponding to the O–O–O bonding intermolecular motion, and the other at ~250 cm⁻¹, associated with the O–O intermolecular stretching. The effect of salt concentration on the vibrational spectrum of water is small and only noticeable in the 1.0 mol·kg⁻¹ CaCl₂ solution, where the two H-bonding modes are shifted to higher frequencies. On the other hand, nanoconfinement substantially modifies the VDOS of water. The intermolecular O–O–O bending and O–O stretching peaks are both blue-shifted, which could be related to changes in the distribution of HBs from a certain tetrahedral order in pure liquid water to a network containing a higher population of two and three HBs per water molecule (Table 2). In Figure 5b, the intermolecular librational (rotational) modes—which occur in the 300–1000 cm⁻¹ region—also display a pronounced spectroscopic shift to higher frequencies as the level of water confinement increases.

Librational modes can be described in terms of motion about the three principal axes of water: wagging, twisting, and rocking. These modes are highly dependent on the local environment and ordering of water [69,70] and sharp librational features have been observed for crystal water in hydrated solids [71]. Therefore, the variations in librational mode frequencies in Figure 5b can be attributed to localized H-bonding and steric restrictions caused by nanoscale confinement.

3.2. Kinetics of Water Exchange at the HAP-Solution Interface

Table 4. Number of accounted water exchange events (N_ex) in the coordination shell of the calcium ions with a duration of more than 0.5 ps. For the HAP-solution systems, the values of N_ex have been normalized with respect to the number of calcium atoms at the surface. Concentration (b) in mol·kg⁻¹.

	b	tsim (ps)	Nex	Nex/100 ps	log kw (s−1 ) a
Solution
Ca2+	–	8000	1532	19.2	11.3
CaCl2	0.6	8000	396.8	5.0	10.7
	1.3	8000	225.1	2.8	10.4
CaF2	0.6	8000	308.1	3.9	10.6
	1.3	8000	202.5	2.5	10.4
HAP, H = 110 Å
Water	–	2000	317.3	15.9	11.2
HAP, H = 60 Å b
Water	–	4000	337.6	8.4	10.9
CaCl2	0.2	4000	374.4	9.4	11.0
	0.5	4000	312.4	7.8	10.9
	1.0	4000	248.6	6.2	10.8
CaF2	0.2	3000	336.9	11.2	11.1
	0.5	3000	296.7	9.9	11.0
	1.0	3000	278.4	9.3	11.0
HAP, H = 40 Å
Water	–	4000	367.4	9.2	11.0
HAP, H = 20 Å
Water	–	4000	111.9	2.8	10.4

^a Water exchange frequency, i.e., the number of exchange events per second. ^b For the HAP nanopores in contact with CaCl₂ and CaF₂ solutions, only water exchanges around surface calcium ions were taken into account.

reports the number of water exchange events (N_ex) in the first coordination shell of the calcium ions, together with the number of water exchange frequency (log k_w)—i.e., the number of exchange events per second—which represents a measure of the ‘lability’ of the hydration shell of the metal [56,58]. Compared with the hydrated calcium ion, the presence of other ions in solution stabilizes the first hydration shell of Ca²⁺ and the number of water exchanges, normalized to 100 ps, reduces from 19.2 in pure liquid water to 2.5–5 exchanges in CaCl₂ (aq) and CaF₂ (aq), depending on the concentration of the salt. The dynamics of water at the HAP–water interface decreases with the degree of confinement. For example, in the HAP nanopore with pore size H = 20 Å, the water exchange frequency (log k_w = 10.4 s⁻¹ ) is one order of magnitude lower than for the hydrated calcium (11.3 s⁻¹). The combined effect of electrolytes in solution (CaCl₂ and CaF_2­) and nanoconfinement was investigated for the HAP (H = 60 Å) nanopore, and the results in Table 4 indicate that the water exchange dynamics at the interface decreases as the concentration of electrolytes in solution increases.

Chloride and fluoride ions impact on the reactivity of calcium by directly coordinating to the metal ion [72] but also when they are in the second coordination shell of Ca²⁺ [61]. In Figure 6, the Ca–O(water), Ca–Cl, and Ca–F RDFs in the HAP (H = 60 Å) nanopores show that surface calcium ions are tightly coordinated to water molecules whereas chloride and fluoride ions are mainly located in the second or higher hydration layers. The only exception are the CaF₂ 0.2 mol·kg⁻¹ and 1.0 mol·kg⁻¹ solutions, where the Ca–F RDFs have a maximum at 2.1 Å, which corresponds to the adsorption of fluoride ions on the surface. In fact, calculation of the HAP surface coverage with respect to the halide ions shows that approximately 1.5% of the calcium ions are coordinated fluoride. On the other hand, the coverage with respect to chloride ions is zero.

The speciation of the calcium–chloride and calcium–fluoride species is also presented in

Table 5. Average coordination number of the surface calcium ions with the water molecules, ${\mathrm{CN}}_{\mathrm{avg}}^{{\mathrm{H}}_2\mathrm{O}}$. Hydroxyapatite surface coverage with respect to the halide ions (X = Cl or F). Speciation of calcium halide ion pairs in solution: contact ion pairs (CIP), solvent-shared ion pairs (SSHIP), solvent-separated ion pairs (SSIP), and free ion pairs (FIP). Concentration (b) in mol·kg⁻¹.

HAP, H = 60 Å	b	$\mathit{C}{\mathit{N}}_{\mathit{a}\mathit{v}\mathit{g}}^{{\mathit{H}}_2\mathit{O}}$	Ca-X Surface Coverage (%)	Speciation of MX Pairs
				CIP	SSHIP	SSIP	FIP
CaCl2	0.2	1.9	0.0	11.5	21.5	24.4	42.6
	0.5	1.8	0.0	18.9	29.4	36.7	15.0
	1.0	1.7	0.1	63.2	19.9	15.4	1.5
CaF2	0.2	2.0	1.5	21.6	27.6	23.2	27.6
	0.5	1.8	0.0	28.9	48.7	19.7	2.8
	1.0	1.8	1.3	28.6	60.1	10.2	1.1

. Fluoride ions are in a more hydrated state than chloride ions and, on average, they are closer to the surface and therefore can have a greater influence on the dynamics of water molecules at the HAP-solution interface. In fact, the number of water exchanges is generally slower when the HAP surface is in contact with CaF₂ solutions (see Table 4). The number of water molecules coordinated to calcium ions (${\mathrm{CN}}_{\mathrm{avg}}^{{\mathrm{H}}_2\mathrm{O}}$) are also reported in Table 5. On average, approximately two water molecules are coordinated to each surface calcium ion.

4. Conclusions

Molecular dynamics simulations of hydroxyapatite nanopores of different sizes (20 Å ≤ H ≤ 110 Å) in contact with water and aqueous electrolyte solutions (CaCl₂ and CaF₂) were conducted in order to determine the role of nanoconfinement and ions in solution on the properties of water and the reactivity of hydroxyapatite surfaces. Analysis of the properties of water has been done via the calculation of pair and angular distribution functions, H-bonding statistics, water self-diffusion coefficients, velocity autocorrelation functions, and power spectra of water. The results have shown that the combined effects of solution composition and nanoconfinement significantly slow the dynamics of water molecules compared with bulk liquid and enhance the rigidity of the H-bonding network. Nanoconfinement, in particular, leads to a strengthening of the H-bonding water network.

The reactivity of hydroxyapatite surfaces was characterized in terms of the dynamics of hydration–dehydration around the calcium sites, which is generally accepted to be the rate-limiting step to crystal growth from aqueous solution of ionic crystals. The dynamics of water at the HAP–water interface decreases substantially with the degree of confinement, and in the HAP (H = 20 Å) nanopore the water exchange frequency is one order of magnitude lower than for the hydrated calcium. Electrolytes in solution also decrease the frequency of water exchange around the calcium ions at the surface, an effect that is dependent upon the specific ion-water interaction.

Together, these changes in the properties of water and surfaces impart an overall rigidifying effect on the water network and reduce the reactivity of the hydroxyapatite-solution interface. Since one of the key indicators of the chemical reactivity of mineral surfaces is the kinetics of water-exchange around the metal cation [28,29,30], atomic-scale information regarding the hydration–dehydration dynamics around the cation can be ‘scaled-up’ to predict the effect on the growth from solution of ionic crystals [58]. This work therefore indicates that nanoconfinement of water within HAP nanopores significantly affects the water exchange reactivity, and the consequent growth kinetics, of the hydroxyapatite surface.