Water Adsorption on the Stoichiometric (001) and (010) Surfaces of

Jan 22, 2009 - Università del Piemonte Orientale “A. Avogadro” and INSTM. ... Lorenzo Zamirri , Marta Corno , Albert Rimola , and Piero Ugliengo ...
1 downloads 0 Views 3MB Size
Download PDF
2188

Langmuir 2009, 25, 2188-2198

Water Adsorption on the Stoichiometric (001) and (010) Surfaces of Hydroxyapatite: A Periodic B3LYP Study Marta Corno,† Claudia Busco,‡ Vera Bolis,‡ Sergio Tosoni,† and Piero Ugliengo*,† Dipartimento di Chimica IFM, Facolta` di Scienze MFN, UniVersita` degli Studi di Torino and NIS Nanostructured Interfaces and Surfaces - Centre of Excellence, Via P. Giuria 7, 10125 Torino, Italy, and Dipartimento DiSCAFF, UniVersita` del Piemonte Orientale “A. AVogadro” and INSTM (Consorzio InteruniVersitario Nazionale per la Scienza e Tecnologia dei Materiali) Unita` del Piemonte Orientale, Largo Donegani 2, 28100 NoVara, Italy ReceiVed October 3, 2008. ReVised Manuscript ReceiVed December 1, 2008 H2O adsorption on hexagonal hydroxyapatite (001) and (010) stoichiometric surfaces has been studied at B3LYP level with a localized Gaussian basis set of polarized double-ζ quality using the periodic CRYSTAL06 code. Because four Ca2+ cations are available at both surfaces, the considered H2O coverages span the 1/4 e θ e 5/4 range. The affinity of both HA surfaces for H2O is large: on the (001) surface, H2O adsorbs molecularly (binding energies BE ≈ 80 kJ mol-1 per adsorbed molecule), whereas it dissociates on the (010) surface, giving rise to new surface terminations (CaOwHw and POHw). The highly negative reaction energy for H2O dissociation (between -250 and -320 kJ mol-1 per adsorbed H2O molecule) strongly suggests that the pristine (010) surface “as cut” from the hydroxyapatite bulk cannot survive in aqueous environment. Conversely, on the reacted surface, H2O adsorbs molecularly with BE similar to those computed for the (001) surface. The B3LYP BEs have been contrasted to the experimental water adsorption enthalpies measured by microcalorimetry on polycrystalline hydroxyapatite samples, showing a fairly good agreement and supporting the suggestion that H2O vapor adsorbs on the already reacted (010) crystalline faces. Harmonic B3LYP vibrational features of adsorbed H2O show, when compared to modes of the gas-phase H2O, a hypsochromic shift of the HOH bending mode (〈∆δHOH〉 ) 49 cm-1) and a bathochromic shift of the OH stretching modes larger than 1700 cm-1 (〈∆νOH 〉 ) 427 cm-1), which are both in good agreement with literature experimental data.

1. Introduction Hydroxyapatite (HA, Ca10(PO4)6(OH)2, apatite family) represents a key mineral in the field of biomaterials for bones and tooth, in terms of repair, reconstruction, and replacement of diseased or damaged tissues. Its relevance is due to the fact that HA is the major component of the mineral phase in mammalian bones and tooth enamel, mostly in its carbonated form (HCA).1-4 Thus, the knowledge of HA surface properties and reactivity becomes an essential requirement to the study of HA-based biomaterials. Knowledge of interfacial atomic structures is then needed for the interpretation of the inorganic-organic interface. The role of H2O in favoring or mediating the interaction between HA surfaces and biomolecules in vivo is still an open field of research and presents great challenge.5,6 To this respect, several experimental and theoretical studies have investigated the HA (001) surface, both bare and in interaction with different molecules. The (001) plane is the dominant surface in the thermodynamic morphology,7-9 but because HA platelets are elongated along the c direction in bones * Corresponding author. E-mail: [email protected]. † Universita` degli Studi di Torino and NIS. ‡ Universita` del Piemonte Orientale “A. Avogadro” and INSTM.

(1) Currey, J. D. Proc. Inst. Mech. Eng., Part H: J. Eng. Med. 1998, 212, 399. (2) Fratzl, P.; Gupta, H. S.; Paschalis, E. P.; Roschger, P. J. Mater. Chem. 2004, 14, 2115. (3) Roveri, N.; Palazzo, B. Hydroxyapatite Nanocrystals as Bone Tissue Substitute. In Tissue, Cell and Organ Engineering; Kumar, C. S. S. R., Ed.; Wiley-VCH: Weinheim, 2006; Vol. 9, p 283. (4) Weiner, S.; Wagner, H. D. Annu. ReV. Mater. Sci. 1998, 28, 271. (5) Kasemo, B. Curr. Opin. Solid State Mater. Sci. 1998, 3, 451. (6) Kasemo, B. Surf. Sci. 2002, 500, 656. (7) Mkhonto, D.; de Leeuw, N. H. J. Mater. Chem. 2002, 12, 2633. (8) de Leeuw, N. H.; Rabone, J. A. L. CrystEngComm 2007, 9, 1178. (9) Vallet-Regı´, M.; Gonza´lez-Calbet, J. M. Prog. Solid State Chem. 2004, 32, 1.

and during biomineralization of tooth hard tissues, the {100} form is the one more developed and responsible for the interaction with molecules,10,11 as shown by HRTEM studies.12 H2O is ubiquitous in the physiological and biochemical environment, and it is involved in the biomineralization processes. Because of the exposed Ca2+ ions at the HA surfaces, H2O is expected to interact as a Lewis base. However, the copresence of electron-rich oxygens belonging to the PO4 groups should also induce the formation of rather strong H-bonds with the adsorbed H2O. In the literature, Zahn et al. have studied the interfaces between HA and H2O by means of classical molecular dynamics simulation.13 The (001) and (001j) surfaces of monoclinic hydroxyapatite (obtained as slices of a block of 3 × 3 × 3 unit cells) were both hydrated by a very large number of H2O molecules. Several strongly ordered hydration layers resulted from the simulation, the first of which is characterized by a network of strong electrostatic/hydrogen bonds that reduce the H2O mobility. As a further outcome, the HA OH groups are not directly protonated by H2O.13 Other studies concerning the interaction of HA and H2O were conducted by de Leeuw,14,15 focused on the dissolution of HA in H2O and the substitution of internal OH- by F-. The interest for OH- substitution by F- ion strengthens the resistance of HA toward dissolution, a fact relevant in the tooth enamel implants.15 To address these points, conventional classical MD simulation procedures with threedimensional periodic boundary conditions were applied on a (10) Simmer, J. P.; Fincham, A. G. Crit. ReV. Oral Biol. Med. 1995, 6, 84. (11) Kirkham, J.; Brookes, S. J.; Shore, S. R.; Wood, R. C.; Smith, D. A.; Zhang, J.; Chen, H.; Robinson, C. Curr. Opin. Colloid Interface Sci. 2002, 7, 124. (12) Sato, K.; Kogure, T.; Iwai, H. J. T. J. Am. Ceram. Soc. 2002, 85, 3054. (13) Zahn, D.; Hochrein, O. Phys. Chem. Chem. Phys. 2003, 5, 4004. (14) de Leeuw, N. H. Phys. Chem. Chem. Phys. 2004, 6, 1860. (15) de Leeuw, N. H. J. Phys. Chem. B 2004, 108, 1809.

10.1021/la803253k CCC: $40.75  2009 American Chemical Society Published on Web 01/22/2009

Water Adsorption on Surfaces of Hydroxyapatite

simulation cell consisting of a 2 × 2 × 3 slab of HA crystal (20.6 Å thick). The (001) slab was artificially replicated to satisfy the periodic boundary conditions, and the empty gap of 50 Å from the next slab was filled with a large number of H2O molecules. The simulation showed some incipient HA dissolution, in that surface Ca ions lengthened their bonds with the lattice oxygen ions due to the solvation of the surface by H2O molecules. Using again classical molecular dynamics techniques, Pan et al. have studied the H2O behavior on hydroxyapatite (100) and (001) crystal faces.16,17 The hydration process envisaged a large number of H2O molecules on both surfaces, and the simulation revealed the formation of three highly structured H2O layers on both faces. Only recently, Ma and Ellis18 reported an ab initio density functional study on the initial stages of the hydration process and Zn substitution on hydroxyapatite (001) surfaces using the PW91 functional and the plane waves/pseudo potential method. Twodimensional slabs with different terminations were considered, and higher H2O loadings were simulated (from 1 to 6 H2O per surface unit cell). The same methodology adopted in ref 14 was followed so that the slabs are replicated using a vacuum gap of about 13 Å. Hydration is studied also in case of fluoroapatite surfaces, both with experimental19,20 and with theoretical methods.7 The most recent paper reporting first-principles calculations based on the SIESTA code (PBE functional) with both stoichiometric (001), (010), (101) and nonstoichiometric Ca- and P-rich (010) surfaces is due to Astala and Stott.21 Calculations showed that the (001) surface was the most stable one. The need to explore nonstoichiometric surfaces stems from the recent HRTEM study by Sato et al.12 in which it was argued that the (010) surface is terminated just over the OH channels running along the (010) plane, defining a Ca-rich termination. Astala and Stott showed, however, that the region of stability of the nonstoichiometric (010) surfaces (both Ca- and P-rich) is outside the window of stability defined by bulk HA, Ca(OH)2, and β-tricalcium phosphate, all needed to define the surface energy formation for nonstoichiometric surfaces. H2O interaction was also investigated and found to be molecular for all but the stoichiometric (010) surface in which H2O dissociates at the surface with the formation of new surface OH functionalities. As it results from the brief summary reported above, only very few studies dealt with water adsorption on the HA most important faces by ab initio computations. For this reason and taking into account the HA relevance in biomaterials field, we have recently started modeling the hexagonal HA bulk22 and its (001) surface23 at B3LYP level with a double-ζ polarized Gaussian basis set and the CRYSTAL06 code.24 The present work aims at extending the B3LYP simulation of this material by investigating the adsorption of H2O on the stoichiometric (001) and (010) faces, which are the HA most relevant exposed faces. The cases of nonstoichiometric (010) surfaces will not be addressed here and will be the focus of future work. Structures, energetic, and vibrational features of the adsorbed H2O have been obtained. (16) Pan, H.; Tao, J.; Wu, T.; Tang, R.; Chin, J. Inorg. Chem. 2006, 22, 1392. (17) Pan, H.; Tao, J.; Tang, R. Front. Chem. China 2007, 2, 156. (18) Ma, X.; Ellis, D. E. Biomaterials 2008, 29, 257. (19) Pareek, A.; Torrelles, X.; Angermund, K.; Rius, J.; Magdans, U.; Gies, H. Langmuir 2008, 24, 2459. (20) Park, C.; Fenter, P.; Zhang, Z.; Cheng, L.; Surchio, N. C. Am. Mineral. 2004, 89, 1647. (21) Astala, R.; Stott, M. J. Phys. ReV. B 2008, 78, 075427. (22) Corno, M.; Busco, C.; Civalleri, B.; Ugliengo, P. Phys. Chem. Chem. Phys. 2006, 8, 2464. (23) Corno, M.; Orlando, R.; Civalleri, B.; Ugliengo, P. Eur. J. Mineral. 2007, 19, 757. (24) Dovesi, R.; Saunders, V. R.; Roetti, C.; Orlando, R.; Zicovich-Wilson, C. M.; Pascale, F.; Civalleri, B.; Doll, K.; Harrison, N. M.; Bush, I. J.; D’Arco, P.; Llunell, M. CRYSTAL2006 User’s Manual; University of Turin: Torino, Italy, 2006; http://www.crystal.unito.it.

Langmuir, Vol. 25, No. 4, 2009 2189

The adoption of the hybrid B3LYP functional for treating semiconductor and oxidic systems has been recently reviewed.25 In the past, it has been proved that B3LYP is a better choice as compared to LDA and well-known GGA functionals concerning the calculation of accurate vibrational features of systems in which H-bond is a key factor.26-28 Microcalorimetric measurements were also carried out to evaluate the strength of the HA/ H2O surface interactions. The experimental adsorption enthalpies measured as a function of increasing water coverage were compared to the B3LYP binding energies. A punctual comparison of available literature experimental data, in particular with IR spectroscopic data from refs 29-31, has also been addressed.

2. Computational and Experimental Details 2.1. Computational. All calculations about hydroxyapatite (001) and (010) surfaces, both free and in interaction with H2O molecules, have been performed with the ab initio CRYSTAL06 code,24 running on standard Linux PC. This code implements the Hartree-Fock and Kohn-Sham self-consistent field method for the study of periodic systems,32 and, in the present version, it allows one to perform full geometry optimization (both internal coordinates and cell size),33 as well as to compute the phonon spectrum at Γ point of molecules, slabs, and crystals.34 Geometrical manipulation and visualization have been carried out with the molecular graphics program MOLDRAW.35 2.1.1. Basis Set. The multielectron wave function is described by linear combination of crystalline orbitals (CO), which, in turn, are expanded in terms of Gaussian-type basis sets. In the present work, calcium ions have been described with the Hay-Wadt small core pseudopotential,36,37 in which the first 10 electrons (shells 1 and 2) are represented by the effective core pseudopotential functions, while for the remaining 10 (shell orbital 3 and 4s), the basis set consists of three contracted Gaussian type functions and one GTF for the outer sp shell (basis [HAYWSC]-31G).38 The exponent of the most diffuse shells is Rsp ) 0.5 bohr-2. Phosphorus atom is described with the basis 85-21d1G with the most diffuse shell exponent Rsp ) 0.135 bohr-2 and Rd ) 0.74583 bohr-2 for polarization, respectively. Finally, oxygen and hydrogen are both represented with a 6-31G* basis set with the most diffuse shell exponents Rsp ) 0.2742 bohr-2 and Rd ) 0.538 bohr-2 for polarization; Rsp ) 0.1613 and Rp ) 1.1 bohr-2, respectively. H2O oxygen and hydrogens were described with the same Gaussian functions used for HA surfaces. 2.1.2. Hamiltonian and Phonon Frequencies. For both (001) and (010) HA surfaces in interaction with H2O, the B3LYP Hamiltonian has been adopted, Becke’s three-parameter (B3) hybrid exchange (25) Cora`, F.; Alfredsson, M.; Mallia, G.; Middlemiss, D. S.; Mackrodt, W. C.; Dovesi, R.; Orlando, R. Struct. Bonding (Berlin) 2004, 113, 171. (26) Merawa, M.; Civalleri, B.; Ugliengo, P.; Noel, Y.; Lichanot, A. J. Chem. Phys. 2003, 119, 1045. (27) Merawa, M.; Labeguerie, P.; Ugliengo, P.; Doll, K.; Dovesi, R. Chem. Phys. Lett. 2004, 387, 453. (28) Ugliengo, P.; Pascale, F.; Me´rawa, M.; Labe´guerie, P.; Tosoni, S.; Dovesi, R. J. Phys. Chem. B 2004, 108, 13632. (29) Bertinetti, L.; Tampieri, A.; Landi, E.; Martra, G.; Coluccia, S. J. Eur. Ceram. Soc. 2006, 26, 987. (30) Bertinetti, L.; Tampieri, A.; Landi, E.; Ducati, C.; Midgley, P.; Coluccia, S.; Martra, G. J. Phys. Chem. C 2007, 111, 4027. (31) Bertinetti, L.; Tampieri, A.; Landi, E.; Bolis, V.; Busco, C.; Martra, G. Key Eng. Mater. 2008, 361-363, 87. (32) Pisani, C.; Dovesi, R.; Roetti, C. Hartree-Fock ab-initio treatment of crystalline systems. In Lecture Notes in Chemistry; Springer: Berlin/Heidelberg/ New York, 1988; Vol. 48, p 193. (33) Civalleri, B.; D’Arco, P.; Orlando, R.; Saunders, V. R.; Dovesi, R. Chem. Phys. Lett. 2001, 348, 131. (34) Zicovich-Wilson, C. M.; Pascale, F.; Roetti, C.; Saunders, V. R.; Orlando, R.; Dovesi, R. J. Comput. Chem. 2004, 25, 1873. (35) Ugliengo, P. MOLDRAW: A molecular graphics program to display and manipulate molecular structures, H1 (32-bit) ed.; Torino, Italy, 2005; http:// www.moldraw.unito.it. (36) Hay, P. J.; Wadt, W. R. J. Chem. Phys. 1985, 82, 299. (37) Hay, P. J.; Wadt, W. R. J. Chem. Phys. 1985, 82, 284. (38) Habas, M. P.; Dovesi, R.; Lichanot, A. J. Phys.: Condens. Matter 1998, 10, 6897.

2190 Langmuir, Vol. 25, No. 4, 2009 functional39 in combination with the gradient-corrected correlation functional of Lee, Yang, and Parr.40 The DFT implementation is such that the exchange-correlation contribution results from a numerical integration of the electron density and its gradient, performed over a grid of points. In CRYSTAL, the Gauss-Legendre quadrature and Lebedev schemes are used to generate, respectively, angular and radial points of the grid. A good compromise between accuracy and cost of calculation for geometry optimization and vibrational frequencies has been shown by using a pruned grid consisting of 75 radial points and 1 subinterval with 974 angular points (75, 974).41 For the largest case treated here, the grid was about 2 800 000 points, which ensured an accuracy on the total number of electrons of 2 × 10-4 electrons on the 900 present in the unit cell. Default values of the tolerances that control the Coulomb and exchange series have been adopted, so that when the overlap between two atomic orbitals is smaller than 10-6 (Coulomb) or 10-12 (exchange), the integral is disregarded. The smaller are these tolerances, the more accurate will be the evaluation of the total energy, because a larger number of bielectronic integrals will be rigorously evaluated at the price of an increasing cost. For the sake of brevity, no more details are provided here, and the interested reader can refer either to the CRYSTAL06 manual24 or to a recent review about the program features.42 The Hamiltonian matrix has been diagonalized43 in 10 reciprocal lattice points (k-points), corresponding to shrinking factor 4.24 For the largest case treated here, the size of the Hamiltonian matrix envisaged 1364 × 1364 elements. Total energy was considered to be converged when the difference between the energy of two subsequent SCF cycles was less than 10-7 hartree. A full relaxation of both lattice parameters and atomic coordinates was performed during the optimization process. The geometry optimization was performed by means of a quasi-Newton algorithm in which the quadratic step (BFGS Hessian updating scheme) is combined with a linear one (parabolic fit) as proposed by Schlegel.44 Convergence was tested on the rms and the absolute value of the largest component of the gradients and the estimated displacements. The threshold for the maximum force, the rms force, the maximum atomic displacement, and the rms atomic displacement on all atoms were set to 0.00045, 0.00030, 0.00180, and 0.00120 au, respectively. The optimization was considered complete when the four conditions were simultaneously satisfied. For a full account of the optimization algorithm, the interested reader can refer to ref 33. Phonon frequencies calculations of H2O adsorbed on HA surfaces have been restricted to molecular fragments, that is, the H2O molecules adsorbed in the unit cell. This is fully justified by the large separation between the water vibrational modes and those of the HA surfaces. Frequencies have been calculated as the eigenvalues obtained by diagonalizing the weighted Hessian matrix associated with the Γ point (point k ) 0 in the first Brillouin zone, called central zone, IR and Raman spectra refer to this point). Second derivatives of the total energy have been computed by finite differences of the analytical first derivative adopting a value of ui ) 0.001 Å as displacement for each atomic coordinate and setting the tolerance of the convergence of the SCF cycles to 10-11 hartree. For more technical points about the implementation and numerical accuracy of the vibrational frequency calculation as implemented in CRYSTAL06, the reader can refer to ref 34. 2.1.3. Binding Energy. In a periodic treatment of surface adsorption processes, the binding energy (BE) per water molecule per unit cell can be defined as the net energy required to remove the adsorbate (here H2O) from the surface to its gas phase, as in the following expression: (39) Becke, A. D. J. Chem. Phys. 1993, 98, 5648. (40) Lee, C.; Yang, W.; Parr, R. G. Phys. ReV. B 1988, 37, 785. (41) Pascale, F.; Tosoni, S.; Zicovich-Wilson, C. M.; Ugliengo, P.; Orlando, R.; Dovesi, R. Chem. Phys. Lett. 2004, 396, 4. (42) Dovesi, R.; Civalleri, B.; Orlando, R.; Roetti, C.; Saunders, V. R. Ab Initio Quantum Simulation in Solid State Chemistry. In ReView in Computational Chemistry; Lipkowitz, K. B., Larter, R., Cundari, T. R., Eds.; John Wiley & Sons Inc.: New York, 2005; Vol. 21, p 1.

Corno et al.

BE ) E(S//S) + EM(W//W) -E(SW//SW)

(1)

where E(S//S) is the energy of the bare slab S in its optimized geometry, EM(W//W) is the molecular energy of the free H2O molecule in its optimized geometry, and E(SW//SW) is the energy of the slab/H2O system in its optimized geometry (the symbol following the double slash identifies the geometry at which the energy has been computed). BE is a positive quantity for a bound adsorbate. For the present case, the deformation energy due to the changes in geometry and the lateral interactions between H2O molecules are easily defined, according to the following equations:

BE ) BE* - δES-δEW ) BE* - δES-∆EW-∆EL (2) δES ) E(S//SW) -E(S//S)

(3)

δEW ) E(W//SW) -EM(W//W) ) ∆EW + ∆EL

(4)

BE* ) E(S//SW) + E(W//SW) -E(SW//SW)

(5)

The different contributions taken into account are (i) δES, the deformation energy of the HA slab, obtained as the difference between E(S//SW), the energy of the bare slab S in the optimized geometry of the whole slab/H2O system SW, and E(S//S) (vide supra); and (ii) δEW, which includes both the deformation energy of the H2O molecule, ∆EW, and the lateral interactions, ∆EL. The first is defined as the purely H2O deformation energy:

∆EW)EM(W//SW) -EM(W//W)

(6)

in which EM(W//SW) is the energy of molecular H2O in the adsorbed geometry, so that ∆EW > 0. The second contribution is due to lateral interactions between the infinite H2O images in the same spatial configuration occurring in the SW periodic system, according to

∆EL)E(W//SW) -EM(W//SW)

(7)

∆EL can be either positive (repulsive) or negative (attractive). With these positions, the BE* ) BE + δES + ∆EW + ∆EL is the binding energy free of deformation and lateral interactions terms. Moreover, dealing with a basis set of localized Gaussian functions, the BE value has to be corrected for the basis set superposition error (BSSE), which causes an overestimation of the BE, when not properly accounted for. Its operational definition is:45

BEC)BE*C-δES-δEW)BE*C-δES-∆EW-∆ELC (8) ∆ELC ) E(N · W//SW)-

∑ j ∑ i(j*i) EM([Wi]Wj//SW)

(9)

BE*C ) (E(S[W]//SW) + E([S]W//SW)) -E(SW//SW) (10) BSSE ) BE - BEC

(11)

The E(S[W]//SW) and E([S]W//SW) terms are the total energy of the slab plus the ghost functions of the water layer and vice versa. The BSSE correction to the ∆EL is important when more than one water molecules is present in the unit cell (N H2O molecules are assumed to be present to define the ∆ELC expression). If not properly taken into account, the BSSE may artificially increment the BE for higher H2O loading cases. The same treatment discussed here has already been used by some of us for the study of the interaction between glycine46 and water47 with silica surfaces. 2.2. Experimental. 2.2.1. Adsorption Microcalorimetry of H2O Vapor. The enthalpy change associated with the adsorption of H2O vapor on a high-surface area hydroxyapatite (HA) specimen was (43) Monkhorst, H. J.; Pack, J. D. Phys. ReV. B 1976, 8, 5188. (44) Schlegel, H. B. J. Comput. Chem. 1982, 3, 214. (45) Boys, S. F.; Bernardi, F. Mol. Phys. 1970, 19, 553. (46) Rimola, A.; Sodupe, M.; Tosoni, S.; Civalleri, B.; Ugliengo, P. Langmuir 2006, 22, 6593. (47) Tosoni, S.; Civalleri, B.; Pascale, F.; Ugliengo, P. J. Phys.: Conf. Ser. 2008, 117, 012026.

Water Adsorption on Surfaces of Hydroxyapatite measured at 303 K by means of a heat-flow microcalorimeter (Calvet C80, Setaram, France) connected to a high vacuum gas volumetric glass apparatus (residual pressure p e 10-5 Torr, 1 Torr ) 133.3 Pa). A well-established stepwise procedure was followed,48-50 which allows one to determine (during the same experiment and for subsequent small increments of the adsorptive) both adsorbed amounts (nads, mol) and integral heats evolved (Qint, J) as a function of the increasing equilibrium pressure (peq). Thanks to the differential construction of the apparatus, all parasite effects (i.e., all effects other than the one due to the interaction(s) of the vapor with the solid surface) were compensated. The adsorptive pressure was monitored by means of a transducer gauge (Ceramicell 0-100 Torr, Varian). Volumetric (nads vs peq) and calorimetric (Qint vs peq) isotherms were obtained,51 but they will not be reported for the sake of brevity, in that they are not useful for the understanding of the present data. Similarly, the integral molar heats of adsorption ([q ) Qint/nads]p, kJ/mol) will not be discussed in that they are comprehensive of the whole thermal contributions arising from the variety of interactions the molecular probe experiences at that peq, and are not appropriate for a comparison with molecular quantities as the computed binding energies (BE). Conversely, by making the derivative of the function, which describes the evolution of the integral heat upon the increasing coverage, the differential heats of adsorption are generated (qdiff ) -∆adsH, kJ/mol), which are still intrinsically average quantities but quantify with a better accuracy the energy of interaction of the molecular probe with the individual sites, as described in ref 48. Such (-∆adsH) quantities, of experimental origin, can thus be conveniently compared to the interaction energy (BEC) of one water molecule with individual model sites, as obtained through the present B3LYP calculations. The measured (-∆adsH) values will be reported in the present work as a function of the increasing coverage, in that the shape of (-∆adsH vs nads) plots depends on, and actually describes, the surface heterogeneity. The extrapolation of the heats plots to vanishing coverage [q0 ) -(∆adsH)0] yields the enthalpy changes associated with uptake on the most energetic sites, expected to be active in the earliest stages of the adsorption process. The nanosized HA sample (BET surface area 80 m2/g) was synthesized in aqueous medium by dropping a solution of H3PO4 in a Ca(OH)2 suspension. The precipitate was then washed, filtered, and dried in air at room temperature as described in refs 30 and 31. Prior to the adsorption experiments, the sample was outgassed overnight at T ) 573 K (p e 10-5 Torr). Such a pretreatment ensured the highest surface density of water-free Ca2+ cations, which were generated by the highest dehydration temperature compatible with the stability of the structure.30

3. Results and Discussion 3.1. Hydroxyapatite (001) and (010) Surface Models. The HA (001) and (010) surfaces were simulated at B3LYP level within the slab approach by selective cuts of the optimized bulk structure described in a previous paper.22 The optimized HA bulk structure has a hexagonal unit cell of 44 atoms [Ca10(PO4)6(OH)2] within the space group P63. The surfaces considered are both perfect planes, without defects (vacancies, steps, or kink sites). For the (001) case, a double layer consisting of 40 atomic layers, 14 Å thick, has been fully optimized within the P3 layer group.23 As for the (010) case, a slab 13 Å thick was fully characterized within P1 symmetry. The (001) surface derived from the hexagonal HA (Figure 1a, parts I and II) exhibits a dipolar character due to the ferroelectric orientation of the OH groups, which is highlighted in Figure 1a, part I. In a previous work,23 the stability of the (001) surface as a function of the thickness has been carefully checked. The structural and electronic characterization of these models, in (48) Bolis, V.; Busco, C.; Ugliengo, P. J. Phys Chem. B 2006, 110, 14849. (49) Bolis, V.; Cavenago, A.; Fubini, B. Langmuir 1997, 13, 895. (50) Fubini, B.; Bolis, V.; Bailes, M.; Stone, F. S. Solid State Ionics 1989, 32-33, 258. (51) Bolis, V.; Busco, C.; Martra, G.; Bertinetti, L.; Corno, M., in preparation.

Langmuir, Vol. 25, No. 4, 2009 2191

Figure 1. Hexagonal hydroxyapatite surfaces: a (001) and b (010). Part I: Views along the lattice vector a. Part II: Unit cell views parallel the ab plane. Atomic labels adopted for the classification of different Ca sites. z axis defines the slab thickness. Table 1. B3LYP Lattice Parameters (Å), Slab Thickness (Å), Unit Cell Surface Area A (Å2), Surface Energy Esurf (J m-2), and Average Geometrical Distances (Å) for the Hexagonal (001) and (010) HA Surfacesa B3LYP

a

b

(001) 9.307 9.307 (010) 6.992 9.261

thickness 13.5 13.0

A

Esurf 〈P-O〉 〈Ca-O〉 〈O-H〉

75.0 1.043 1.551 64.7 1.709 1.551

2.385 2.370

0.9705 0.9692

a The slab thickness is defined as the perpendicular distance between the most exposed Ca2+ surface ions.

terms of surface energy, geometry relaxation, band gap, field across the slab, and electrostatic features in close proximity of the surface, has shown that the OH ferroelectricity does not prevent the formation of a (001) slab of thickness of at least 10 nm. The HA monoclinic phase is a nonferroelectric alternative to the hexagonal one, but the computational cost rises dramatically because its unit cell contains twice as many atoms as the hexagonal HA. For these reasons, the hexagonal HA double layer slab is considered here as a proper model of the HA (001) face to study H2O adsorption processes. Table 1 reports some significant features of the (001) double layer model, including its surface energy Esurf ) 1.043 J m-2, defined as

Esurf ) (EsN-N · Eb)/2A

(12)

in which EsN, N, Eb, and A are the energy of the slab unit cell, the number of bulk unit cell present in each considered slab model, the energy of the bulk, and the unit cell area, respectively.

2192 Langmuir, Vol. 25, No. 4, 2009

In the literature, different values of the (001) Esurf are reported, depending on the simulated model and on the computational approach. Rulis et al.52 calculated 0.871 J m-2 for a HA (001) supercell slab model (VASP and OLCAO codes, GGA PBE, Ecutoff ) 600 eV). Filgueiras et al.53 have used the METADISE code to run energy minimizations of two HA (001) surfaces, the first terminating as PO4-Ca-PO4-PO4-Ca-Ca with an Esurf of 1.08 J m-2 and the second terminating as Ca-Ca with an Esurf of 1.12 J m-2. de Leeuw and Rabone8 have computed for a simulation cell containing a 2 × 2 × 3 slab of hydroxyapatite a lower value of Esurf, 0.80 J m-2, due to the inclusion of the temperature (310 K) in the classical MD simulations with the DL_POLY code. The other HA surface considered here is the stoichiometric (010) face (see Figure 1b, parts I and II), as it is very extended in the crystal morphology.54 The considered model envisages 88 atoms in the unit cell with a 13 Å slab thickness, so that it is comparable in size to the (001) double layer. The B3LYP (010) surface energy is 1.709 J/m2, higher than that of the (001) double layer (see Table 1), because the number and type of bonds that should be cut to define the (010) surface impart a higher cost in comparison with the (001) surface. Therefore, the (010) model is expected to be more active toward the adsorption of H2O molecules than the (001) one. Optimized lattice parameters, thickness, area, Esurf, and the most significant averaged bond lengths of the (010) model are reported in Table 1. Contrary to what happened for the (001) double layer, for the (010) case the OH dipoles are aligned along the a axis (perpendicularly to the z coordinate, which defines the slab thickness) so that the surface is nonpolar. This is also shown by the electrostatic potential map (not shown) computed above and below the (010) slab, which does not exhibit any macroscopic dipole across the slab. Computed band gap is 6.4 eV, to be compared to the value of 8.1 eV for the (001) model. Our Esurf values are in good agreement with those computed by Astala and Stott21 of 1.20 and 1.68 J/m2 with PBE and the SIESTA code. For both considered models, the analysis of the electrostatic potential features was useful to predict the possible adsorption sites. Indeed, from electrostatic potential maps taken at 2 Å from the top/bottom surfaces (not shown for brevity), highly positive values of potential on top of the Ca2+ ions were located. The same maps also revealed deep negative zones in the proximity of the oxygen atoms of the PO4 groups. The docking of water toward the surface will then follow the “electrostatic complementarity principle”: on Ca2+ ions, H2O behaves as a Lewis acid interacting through the oxygen atom, whereas it H-bonds with the oxygen atoms of the PO4 groups. Because H2O is polarized by the interaction with the Ca2+ ions, the H-bonds’ strength will, in turn, be enhanced as compared to that of free water. The role of Ca2+ ions as polarizing centers was also suggested by Rulis et al.52 who concluded that the (001) surfaces of HA are mostly positively charged due to the presence of calcium ions so that they can be implicated in the adsorption of H2O and other organic molecules in an aqueous environment, relevant for HA bioactivity. In the following, the cases of the two surfaces interacting with H2O will be described separately to improve the readability of the results. For each model, H2O has been adsorbed on both exposed faces (upper and lower), to avoid symmetry or electrostatic artifacts. Once H2O molecules were docked to the (52) Rulis, P.; Yao, H.; Ouyang, L.; Ching, W. Y. Phys. ReV. B 2007, 76, 245410(15). (53) Filgueiras, M. R. T.; Mkhonto, D.; de Leeuw, N. H. J. Cryst. Growth 2006, 294, 60. (54) Wierzbicki, A.; Cheung, H. S. J. Mol. Struct. (THEOCHEM) 2000, 529, 73.

Corno et al.

surface as described above, full B3LYP geometry optimization was carried out to locate minimum structures. In the following, each structure is characterized by means of its geometrical features, average BE with respect to H2O adsorption, and harmonic vibrational frequencies limited to the modes of the adsorbed H2O. 3.2. H2O Adsorption on HA (001) Surface Model. 3.2.1. Geometrical Features of H2O Adsorbed on Different Ca2+ Sites at Increasing CoVerage. Considering the upper face as reference, four Ca ions are available to H2O adsorption (Ca1, Ca3A, Ca3B, and Ca3C) in the HA (001) unit cell (see Figure 1a, parts I and II). The same obviously applies to the lower face. An apex has been used to distinguish between sites at the top and bottom faces; that is, Ca1′ refers to calcium at the bottom face and so on. Ca2 has been excluded as an adsorptive site because it is almost buried in the structure and is considered only as a geometrical reference point inside the cell, to facilitate the description of higher H2O loading. In the literature, calcium ions at the (001) surface are also labeled as CaI/Ca1 and CaII/Ca2, by referring to HA P63/m bulk crystal.17,18,52,55,56 In the present work, CaII/Ca2 should read Ca3, while CaI/Ca1 includes both Ca1 and Ca2, due to the symmetry lowering of the bulk structure when passing from P63/m to P63 space groups (for a deeper explanation, see ref 22). Three different coverages θ ) H2O/ Ca2+surf have been studied considering Ca2+surf ) 4 for each face (see Figure 1, part I, and discussion above), θ ) 0.25 (001_W2/ Ca1 and Ca3, Figure 2), θ ) 0.5 (001_W4, Figure 3, part II), θ ) 1 (001_W8, Figure 3, part III), and θ ) 1.25 (001_W10, Figure 3, part IV). The lowest coverage case (indicated from now on as 001_W2) envisages the adsorption of two H2O molecules per unit cell, one per each face, respectively, on Ca1 (the most exposed cation in the bare surface, its position being very close to that of the “as-cut” surface from HA bulk) and Ca3 (less exposed, forming a triangle surrounding the OH group). H2O molecules were placed with their oxygen Ow at about 2 Å from the chosen calcium ions, Ca1, Ca1′ or Ca3A, Ca3′A as in Figure 2. When only one H2O molecule is present, it easily distinguishes between the two differently exposed cations (four different situations including Ca1′ and Ca3′A). Each H2O molecule interacts with the Ca2+ surface ions through its oxygen (the shortest Ca · · · Ow being for Ca1, 2.293 Å, the longest for Ca3′A, 2.444 Å) and by one H-bonding between the Hw and the surface oxygen Os of superficial phosphate groups. H-bondings are in general quite short, the shortest being 1.51 Å between Hw and Os for the case of adsorption on Ca1 (see Figure 2 for further details). Cases envisaging Ca3B and Ca3C are not reported for the sake of brevity, exhibiting almost the same features of the Ca3A case. Figure 3 shows four models of 001_W adsorption at increasing H2O coverage, viewed perpendicularly to the ab unit cell plane for both upper and lower faces of the surface. The 001_W2 case for Ca1 and Ca1′ (already shown in Figure 2) is repeated to facilitate the comparison. As a whole, the maximum H2O loading is reached with 10 adsorbed H2O, that is, one on each of the four available calcium ions per cell (1 Ca1 + 3 Ca3) at both upper and lower faces of the slab and the extra two H2O molecules adsorbed in the region identified by the position of Ca2 and Ca2′, at a Ca-Ow distance of 4 Å. The computed maximum amount of H2O compares well with the experimental value of 4.5 ( 0.5 H2O molecules per nm2 reported by Bertinetti et al. on the basis of experimental measurements,30 which in turn is in agreement (55) Pan, H.; Tao, J.; Xu, X.; Tang, R. Langmuir 2007, 23, 8972. (56) Chappell, H.; Duer, M.; Groom, N.; Pickard, C.; Bristowe, P. Phys. Chem. Chem. Phys. 2008, 10, 600.

Water Adsorption on Surfaces of Hydroxyapatite

Figure 2. Unit cell views along the b axis for the two 001_W2 models with 2 H2O molecules per unit cell adsorbed on Ca1 or Ca3A (upper and lower surfaces). I: Side view of the unit cell; local view of the adsorption site at upper (II) and lower (III) surfaces. Ow and Os as H2O and surface oxygen, respectively. H-bonds as dotted black lines. Ca ions in cyan, to be distinguished from H2O hydrogen atoms. z axis defines the slab thickness.

with the datum of 4.3 Ca2+ exposed at the surface of hydroxyapatite particles per nm2 measured by Kukura et al.57 The W4 model envisages two H2O molecules adsorbed to each face to give θ ) 0.5 (001_W4, Figure 3, part II). Ca1 and Ca3C coordinate the oxygen atoms of the two water molecules, which are also in mutual H-bond contact. Interestingly, each H2O molecule (either at the top or at the bottom face) establishes rather short H-bond contacts with the basic oxygens belonging to the PO4 anions. For the W8 model (H2O on Ca1 and on the 3 Ca3), four H2O are adsorbed on each face reaching a θ ) 1 coverage. The optimized structure reveals a rather complex network of H-bond interactions between the adsorbed H2O molecules, which increases the lateral interaction contribution to the BE (vide infra). Indeed, H2O molecules tend to cluster via H-bonding inside the unit cell, keeping, at the same time, the Ow as close as possible to the most exposed Ca2+ ions (see Figure 3, parts III and IV for details). Present results also show that the structural exposed OH groups (57) Kukura, M.; Bell, L. C.; Posner, A. M.; Quirk, J. P. J. Phys. Chem. 1972, 76, 900.

Langmuir, Vol. 25, No. 4, 2009 2193

Figure 3. Views of the unit cell on top (left side) and at the bottom for (001) HA double layer models in interaction with 2 (W2/Ca1, I), 4 (W4, II), 8 (W8, III), and 10 (W10, IV) H2O molecules, respectively. Unit cell borders in blue, H2O atoms rendered as balls, Ca ions in cyan, H-bond lengths in angstroms.

of the HA (001) surface are not involved as adsorption sites: H2O moves away from the surface hydroxyl, in agreement with similar findings by Ma and Ellis.18 Only at high θ (001_W8 and 001_W10 models, Figure 3, parts III and IV), one H2O is oriented with its Ow toward the H of the superficial OH group, at a distance of about 2.00 Å, so that the HA O-H distance changes slightly from 0.964 to 0.974 Å. Detailed analysis of the geometrical features of the 001_W8 model (see Figure 3, part III) shows an average Ca-Ow of 2.411 Å (shortest and longest at 2.312 (Ca1) and 2.463 Å (Ca3B′)), to be compared to the corresponding values of 2.586 Å (shortest 2.304 (Ca1) and the longest 3.865 Å (Ca3B)) for the 001_W10 model (Figure 3, part IV). The first adsorption layer resulting from the simulation by Ma and Ellis18 contains 5 H2O located at about 1.6 Å from the surface (considering as the zero the height of oxygen ions belonging to the surface phosphate groups), while the sixth H2O is floating at 3.2 Å from the surface plane. For our highest H2O coverage 001_W10 model (Figure 3IV), all five adsorbed H2O molecules are at a maximum height of about 1.7 Å from the surface plane. The most remarkable result of the present simulation is that for all considered models, H2O adsorbs on HA (001) molecularly and never dissociates, in agreement with similar findings by Ma and Ellis18 and by Astala and Stott.21 Further checks about this point were also carried out, in which H2O was manually placed on the HA (001) surface in a predissociated state by forming Ca-OH and P-O-H groups. In all cases, few optimization

2194 Langmuir, Vol. 25, No. 4, 2009

Corno et al.

Table 2. B3LYP Binding Energies (kJ mol-1 per adsorbed molecule) Corrected for the BSSE, BEC, and Its Different Energy Contributions as a Function of the Increasing Amount of Adsorbed H2O Molecules per Unit Cell (from W2 to W10) for the (001) Slab Modelsa (001) model W2 (Ca1) W2 (Ca3A) W4 W8 W10 a

δES δEW ∆EW ∆EL

∆ELC BE* BE*C BE

BEC

29.9

7.4 8.1

-0.7 -0.7 185.4 138.4 148.1 100.4

17.0

0.3 0.9

-0.6 -0.6 111.6 58.5 94.3 40.6

22.5 1.3 9.0 -7.7 -3.4 151.3 103.5 127.5 75.4 12.7 -15.9 5.0 -20.9 -13.0 117.3 73.3 120.5 68.6 11.2 -25.0 1.5 -26.6 -15.9 102.6 63.4 116.4 66.5

Note that ∆EL ) δEW - ∆EW.

cycles were enough to restore the molecular integrity of the H2O molecule. 3.2.2. Interaction Strength: Binding Energy. Table 2 reports the binding energy BE values together with the BSSE corrected values BEC and their different contributions for the five 001_W models described so far (see computational details section). BE is averaged over the number of H2O molecules adsorbed per unit cell for each considered θ. BSSE turned out to be sizable as can be computed by subtracting the last two columns of Table 2 (see eq 11). Previous experience with basis set quality46,47 similar to that adopted here showed that BEC are reliable because they represent much better than uncorrected BE values the limiting BE value, which would be computed with better quality basis sets. BEC values are a function of the Ca2+ exposure, and, indeed, the strongest interaction refers to the most exposed Ca1 ion (W2/ Ca1Ca1′ case, about 100 kJ mol-1) as shown in Figure 1a for the bare surface23 and in Figure 2 for the adsorbed water case. As expected, much lower BEC values (about 40 kJ mol-1) have been computed on less exposed Ca3 and Ca3′ ions. BEC changes also with H2O loading, in that a decrease with the increasing H2O loading is computed when passing from 2 to 10 adsorbed H2O molecules (BEC from 100.4 (Ca1) to 66.5 for W10). As a global BEC average, a value of 78 kJ mol-1 per molecule resulted at B3LYP, while Ma and Ellis computed 68 kJ mol-1 per molecule at the PBE level.18 In their case, a lower value of 75.3 kJ mol-1 per molecule was computed for the case of one H2O per unit cell, which decreases to 64.5 kJ mol-1 per molecule when six H2O were adsorbed. Table 2 shows the partitioning of the BE in terms of the energy contributions discussed in the computational method section. Interestingly, the deformation energy of HA surface, δES, decreases from 30 to 11 kJ mol-1 along the increase in θ, because water interacts less strongly with the HA surface at high θ and prefers to mutually H-bond, increasing the importance of ∆ELC, which becomes negative (attractive) and large (about -16 kJ mol-1 per molecule) for the highest θ. The H2O deformation term ∆EW is large for both W2/Ca1 and W4 cases, in which a strongly localized interaction occurs with the most exposed Ca ions, as shown by the high δES term. For the W2/Ca3A case, both ∆EW and δEW become small as a consequence of a longer Ca---Ow distance (see previous section) when compared to the W2/Ca1 case. Accordingly, the δES term is about one-half of the value resulting for the W2/Ca1 case. 3.3. H2O Adsorption on HA (010) Surface Model. The same approach described so far for the (001) HA surfaces has been followed to study H2O adsorption on the (010) HA surface, so that in the following the results will be discussed in a similar manner.

Figure 4. Top (I) and local Ca1 adsorption site (II) views for 010_W2/ Ca1. Ca ions in cyan.

3.3.1. Geometrical Features of Adsorbed H2O on Different Ca2+ Sites at Increasing CoVerage. For the HA (010) slab, four different Ca2+ ion types are available for adsorption, that is, Ca1, Ca2A, Ca2B (almost equivalent for symmetry), and Ca3 on the upper face and the equivalent ions labeled with the apex for the lower face (see Figure 1b, parts I and II). In principle, four different coverages are possible, that is, 0.25 e θ e 1. Starting with the θ ) 0.25 case, four 010_W2 models have been optimized at B3LYP level, W2/Ca1Ca1′, W2/Ca2Ca2′A, W2/Ca2Ca2′B, and W2/Ca3Ca3′. For reasons of brevity and because for the (010) surface the upper and lower faces are geometrical and electrically identical, from now on only Ca of the upper surface will be reported. The first dramatic difference for the adsorption on HA (010) face as compared to the (001) case is that H2O spontaneously dissociates on Ca2A, Ca2B, and Ca3, whereas it remains molecularly adsorbed only for the W2/ Ca1 model (see Figure 4). The same finding has been reported by Astala and Stott.21 Dissociation brings about new surface functionalities, that is, CaOwHw2 and POs1Hw1, as shown in Figure 5 for W2/Ca3 (the other two cases are not shown for brevity). Figure 4, parts I and II show the W2/Ca1 structure, in which H2O is adsorbed molecularly, similarly to the (001) face. Some of the most relevant geometrical features are reported in Table 3, included the two H-bonds formed with the Os1 and Os2 surface oxygen of the most exposed phosphate groups. Interestingly, some experimental evidence from the NMR technique that protonation of the fluoroapatite surface may indeed occur (as a function of pH) has been reported.58 The dissociative H2O adsorption is shown in Figure 5 for the W2/Ca3 case and reveals that the newly formed OwHw2 group is shared by Ca3, Ca2A, and Ca2B ions with a local geometry resembling that surrounding the OH group of the HA (001) surface (Figure 5, structure III). It is clear that this particular geometrical arrangement gives stability to the final structure, and it is worth (58) Jarlbring, M.; Sandstrom, D. E.; Antzutkin, O. N.; Forsling, W. Langmuir 2006, 22, 4787.

Water Adsorption on Surfaces of Hydroxyapatite

Langmuir, Vol. 25, No. 4, 2009 2195

Figure 5. B3LYP geometry optimization for 010_W2/Ca3 as seen along the z axis. Starting structure (I), proton sharing situation (II), and final point with dissociated H2O (III). Local view of Ca3 adsorption site at the optimized structure (IV). Labels refer to Table 3 and are in blue for H2O atoms. Dissociative process occurring on lower face not shown for brevity. Table 3. Most Relevant Computed Bond Lengths and Angles for H2O Adsorbed on the Four Different 010_W2 Models, Expressed in angstroms and deg, Respectivelya (010) model W2/Ca1

W2/Ca2A W2/Ca2B W2/Ca3 a

Ca1-Ow Hw1 · · · Os1

Hw2 · · · Os2

Hw1OwHw2 POs1Hw1

2.373

1.863

2.021

101.2

CaOw-Hw2

POs1-Hw1

Os1Hw1 · · · Ow/s2

CaOwHw2

POs1Hw1

0.964 0.964 0.965

1.030 1.029 0.989

1.635 1.693 1.962

123.3 123.5 122.2

109.9 109.9 111.9

Labels refer to Figure 5.

noting that H2O does not dissociate on Ca1 because the surface OH cannot be similarly stabilized. Figure 5 shows three snapshots resulting from the geometry optimization (I starting point, II middle - proton sharing, III final product) for the W2/Ca3 model. Corresponding geometrical data of Table 3 reveal that the newly formed CaOw-Hw2 bonds are very close to the typical CaO-H bond found in the B3LYP optimized structure of Ca(OH)2 portlandite (CaO-H about 0.968 and 0.965 Å for portlandite59 and the present case, respectively). For the (010) surface, higher H2O loadings were also studied to assess if the dissociation of H2O on Ca3 would further promote H2O dissociation on Ca1, which was computed as the only nondissociative site at θ ) 0.25. Model 010_R_W2/Ca1 of Figure 6, part I (the capital R in the model label indicates adsorption on the H2O reacted surface), which corresponds to θ ) 0.25 on (59) Tosoni, S.; Pascale, F.; Ugliengo, P.; Orlando, R.; Saunders, V. R.; Dovesi, R. Mol. Phys. 2005, 103, 2549.

Figure 6. B3LYP unit cell views along the z axis for the 010_R_W2/ Ca1 (I), 010_R_W4/Ca1Ca2A (II), 010_R_W6/Ca1Ca2ACa3 (III), and 010_R_W8/Ca1Ca2ACa2BCa3 (IV) models. H2O labels in blue.

the 010_R surface, shows that H2O on Ca1 remains a molecularly adsorbed spectator for all considered stages of the adsorption process, with no influence on the dissociation occurring on Ca3. Molecularly adsorbed H2O is involved in a rather strong H-bond (OwHw · · · Os ) 1.89 Å) with the oxygen atom of the surface PO4 group (see Figure 6, part I). This is again in nice agreement with the findings by Astala and Stott21 in which the second H2O molecule remains molecularly adsorbed. Three higher coverages were also considered starting from the reacted (010) surface (010_R): (i) two H2O molecules for each face at Ca1 and Ca2A type ions (θ ) 0.5, Figure 6, part II); (ii) three H2O molecules for each face, located on Ca1, Ca2A, and Ca3 ions (θ ) 0.75, Figure 6, part III); and (iii) four H2O for each face, the fourth one located at Ca2B type ions, so reaching the same loading as the 001_W8 model (θ ) 1, Figure 6, part IV). At higher θ, at variance with what was predicted on the (001) face, the lateral interactions between adsorbed H2O are expected to be less important, as direct H-bond contact only involves pairs of H2O (see Figure 6). 3.3.2. Interaction Strength: Binding Energy. The dissociative adsorption of H2O is a barrierless process (geometry spontaneously evolves from structure I toward IV of Figure 5 during optimization) and characterizes a highly spontaneous process, the reaction energy being -315 kJ mol-1 per molecule and -249 kJ mol-1 per molecule for H2O reacted on Ca3 (Figure 5, part III) and Ca2A/Ca2B ions, respectively. The energetics of these dissociative processes has been confirmed by single point energy calculations using a polarized triple-ζ quality basis set.60 Astala and Stott21 computed a value of about 292 kJ/mol using PBE and SIESTA code, in very good agreement with our data derived at (60) Schafer, A.; Horn, H.; Ahlrichs, R. J. Chem. Phys. 1992, 97, 2571.

2196 Langmuir, Vol. 25, No. 4, 2009

Corno et al.

Table 4. Energy Contributions to the Final BEC Value for 010_W2 and 010_R Models, Where H2O Is Molecularly Adsorbed, Expressed in kJ mol-1 per adsorbed moleculea (010) model

δES

δEW

∆EW

∆EL

∆ELC

BE*

BE*C

BE

BEC

W2/Ca1 R_W2/Ca1 R_W4/Ca1Ca2A R_W6/Ca1Ca2ACa3 R_W8/Ca1Ca2A2BCa3

15.4 25.8 25.2 23.1 22.0

1.5 6.9 -2.9 -11.3 -11.4

2.1 5.9 6.4 3.2 2.8

-0.6 0.9 -9.3 -14.4 -14.2

-0.6 0.9 -3.5 -8.0 -4.2

184.4 181.8 163.9 135.6 122.9

134.7 123.9 110.0 88.0 79.6

167.5 149.1 141.6 123.8 112.3

117.8 92.1 81.9 69.7 59.0

a

Note that ∆EL ) δEW - ∆EW.

B3LYP with CRYSTAL06. This shows that, irrespective of the adopted functional or periodic code, the underlying physics remains the same. Table 4 reports the computed BEC and its components for the case of molecular H2O adsorption on Ca1 (Figure 4) and for the other models of molecular adsorption on the reacted (010) surface, that is, after H2O dissociation on Ca3 (010_R/W2, W4, W6, and W8 models of Figure 6, parts I, II, III, and IV, respectively). The resulting BEC values (92, 82, 70, and 59 kJ mol-1 per molecule, respectively) are smaller than that for the molecular adsorption on Ca1 (118 kJ mol-1 per molecule, Figure 4) and decrease with the increase of θ. A BEC of 59 kJ mol-1 per molecule corresponding to θ ) 1, that is, the maximum (010) surface coverage simulated, to be compared to the value of 69 kJ mol-1 per molecule for the (001) crystal face at the same θ, indicates that on the reconstructed surface the affinity for H2O is still high. A detailed analysis of the various contributions to the BE (see Table 4) shows that the geometrical deformations δES suffered by the (010) surface are all higher than those computed for the (001) surface. The lateral interactions, ∆ELC, increase at higher θ, but less than for the (001) case as was already anticipated in the discussion on the geometries (vide supra). 3.4. Vibrational Frequencies of Adsorbed H2O/HA (001) and (010). B3LYP harmonic vibrational frequencies limited to the modes of the adsorbed H2O molecules have been computed on both (001) and (010) faces and compared to those of free H2O. Figure 7 shows the values of the shift of the O-H stretching modes and the HOH bending mode with respect to those computed for the free H2O (3798 and 3918 cm-1 for the OH symmetric and asymmetric stretching modes; 1658 cm-1 for the HOH bending mode). Data in Figure 7 were computed considering that when H2O is adsorbed at the surfaces and H-bonds are formed, it is no longer possible to distinguish between asymmetric and symmetric O-H stretching modes. As an arbitrary choice, modes involved in H-bonds are compared to free H2O symmetric stretching, whereas those free from H-bond are compared to the H2O asymmetric stretching. The general trend consists of a prevalent hypsochromic shift of the HOH bending mode (average value of 49 cm-1), as expected for H-bonded systems61 because the HOH force constant is strengthened due to the hindered motion of the proton engaged in the H-bond. For the OH stretching modes, a large bathochromic shift is predicted in all cases, with an average value of 427 cm-1 and maximum value of 1651 cm-1 for the 001_W8 model. The large OH bathochromic shifts are due to the formation of several H-bonds, as described in the previous section. Indeed, when the bathochromic shift is at its maximum (∼1650 cm-1), water H-bonds with an oxygen of the HA (001) surface, and the corresponding Hw · · · Os bond length is as small as 1.3 Å (model 001_W8, Figure 3, part III, for H2O adsorbed on Ca1, Ow · · · Ca1 ) 2.300 Å, Hw-Ow 1.060 Å). (61) Jeffrey, G. A. An Introduction to Hydrogen Bonding; Oxford University Press: New York, 1997.

Figure 7. IR harmonic frequencies shifts of adsorbed H2O with respect to free H2O for OH stretching symmetric and asymmetric (first two bars) and HOH bending, expressed in cm-1. Light blue, 001_W2/Ca1; black, 001_W4; blue, 001_W8; and red, 001_W10 models; while light green, 010_W2/Ca1; orange, 010_R_W2/Ca1; and pink, 010_R_W4/Ca1Ca2A models.

As already mentioned, for OH stretching modes two cases are present: (i) the O-H stretching in which the H is involved in a H-bond with the surface Os (indicated as symmetric in Figure 7) and (ii) the O-H stretching in which the H is involved in H-bond with nearby H2O molecules (indicated as asymmetric in Figure 7). Remarkably, considering the complexity of the present systems, a correlation between the strength of the H-bond, represented by the OwHw · · · Os distance, and the bathochromic shift of the corresponding stretching mode is shown to hold (see Figure 8). As expected,62 a plateau value in the OH stretching shift is reached at large OwHw · · · Os distances for both (001) and (010) cases (parts a and b of Figure 8, respectively). It is worth noting that all OH stretching shifts are reported, so that also cases without H-bond formation are considered as shown by OwHw · · · Os distances as large as 3.2 and 4 Å, for (001) and (010) surfaces, respectively. Computed B3LYP IR frequencies were tentatively compared to recently published30 experimental data on powdered hydroxyapatite nanocrystals. Experimentally, a broad band extending from 3740 to 2500 cm-1 is seen as a clear indication of H-bond of various strength occurring in correspondence to H2O vibrational stretching modes, due to the different possible sites of adsorption and to different H-bonds strength. From the experiment, the HOH bending mode undergoes a hypsochromic shift, whose average value is about 40 cm-1, in remarkable (62) Libowitsky, E. Monatsh. Chem. 1999, 130, 1047.

Water Adsorption on Surfaces of Hydroxyapatite

Langmuir, Vol. 25, No. 4, 2009 2197

Figure 9. Calorimetrically measured enthalpy adsorption [-∆adsH] reported as a function of water coverage on a nanosized hydroxyapatite specimen (HA, O). Dashed line: water enthalpy of liquefaction [-∆liqH ) 44 kJ mol-1]. Computed binding energies [BEC] values (at proper θ, see the text) are reported for comparison purposes for both (001) and (010) models (red open and blue solid stars, respectively).

Figure 8. Correlation between hydrogen bond lengths OwHw · · · Os and stretching ∆ν(O-H) shifts with respect to free H2O for models involving (001) and (010) faces.

agreement with the average B3LYP value of 49 cm-1. For the OH stretching, average and maximum bathochromic shifts of about 400 and 1100 cm-1 are observed to be compared to corresponding B3LYP values of 427 and 1650 cm-1, respectively. Despite the good general agreement, it has to be considered that the B3LYP values are computed within the harmonic approximation, which gives smaller bathochromic shifts in comparison with those resulting from an anharmonic treatment. However, the presence of a relatively large BSSE in the BE is expected to cause an overestimation of the harmonic OH bathochromic shifts, so that the present values are the result of some error compensation. The CaO-H stretching frequencies associated with the new Ca-OH functionality derived from H2O dissociation (010_W2/ Ca3Ca3′ model) are computed at 3843 (Ca3) and 3874 (Ca3′) cm-1, respectively, in good agreement with the corresponding O-H stretching value of the superficial OH belonging to the (001) face.23 The similarity in the OH stretching values remains so also at higher θ values. The other newly formed surface functionality derived by the H2O dissociation is the POH group. PO-Hw stretching and POHw bending modes are influenced by the formation of H-bond with the oxygen of the new surface functionality CaOwHw. In all models, at different H2O coverages, PO-Hw stretching modes are computed in a range between 3700 and 2700 cm-1, while

POHw bending frequencies span the range 1000-1400 cm-1, corresponding to PO-Hw · · · HwOwCa H-bonding varying between 2.0 and 1.6 Å, respectively. 3.5. H2O Adsorption on Ca2+ Sites: Microcalorimetry versus Computed Binding Energies. In Figure 9, the calorimetrically measured enthalpy of adsorption [-∆adsH] of water vapor at the surface of a nanosized HA specimen is reported as a function of increasing coverage, which in the present case has been expressed as the number of H2O molecules adsorbed per surface cation Ca2+.30 Computed BEC values, as derived from Tables 2 and 4, are reported for comparison purposes for both (001) and (010) model surfaces. Let us first consider the experimental data reported in Figure 9. The enthalpy versus coverage plot shown in the figure is typical of a highly heterogeneous surface.48 The enthalpy of adsorption decreases from a rather high zero-coverage value [(-∆adsH)0 ≈ 120 kJ/mol] down to a value as low as ∼50 kJ/mol. It is worth noting that the enthalpy curve eventually approaches the value of the enthalpy of liquefaction of water [-∆liqH ) 44 kJ mol-1] only after a second water molecule has been adsorbed per Ca2+ cation (as average), suggesting that the H-bonding interactions responsible for the adsorption of water in the second layer are still stronger than those experienced by H2O molecules in the liquid. This effect is most likely due to the strong polarization of the first-layer water adsorbed molecules, which are tightly bound to the coordinatively unsaturated surface cations (cus), as witnessed by the high zero-coverage enthalpy of adsorption. This high enthalpy value is not, however, attributable to a (partial) dissociation of water at the surface sites of our HA sample in that it has been demonstrated by IR spectroscopy (as reported and extensively discussed in ref 30) that, despite their extremely high hydrophilicity, the experimental nanosized surfaces do not dissociate water. The reasons for such a strong affinity for water stem most likely on the synergy between cus Ca2+ cations and PO4 species located in the close vicinity. This is confirmed by the B3LYP optimized structures (Figures 2, 3, and 6) in which strong H-bonds with the PO4 species give an extra stabilization to H2O molecules coordinated on cus Ca2+ sites. To make the proper comparison between measured [-∆adsH] and calculated [B3LYP BEC] energetic data, it is necessary to make a preliminary comment. As already described, the B3LYP results indicate a rather different mechanism for H2O adsorption on the (001) and (010) surfaces: for the former surface, H2O adsorption resulted

2198 Langmuir, Vol. 25, No. 4, 2009

associative (molecular) in nature, whereas on the latter one it was dissociative (more than 250 kJ mol-1 are developed per H2O molecule) and not activated. For these reasons, one can infer that the (010) surface does not survive as such in contact with an aqueous solution, which is typically employed during the synthesis of crystalline HA materials (vide supra experimental section). As a consequence of that, the (010) face depicted in our starting computer model (vide supra Figure 1b) is not really representative of the experimental (010) surface of the HA nanocrystals synthesized in a water-rich environment. In conclusion, for our comparison purposes, instead of the B3LYP BEC values computed for the “unreacted” (010) surface (Figure 1b), we have taken the corresponding values computed for the “reacted” (010) crystalline surface (see Figure 5, part III and data reported in Table 4). This means that the BEC reported in Figure 9 for the (010) surface represent the energy of interaction of H2O molecularly adsorbed at the sites formed upon dissociation and which behave similarly to the (001) surface sites. For both (010) and (001) surfaces, the computed values span in the 100-60 kJ mol-1 range (per adsorbed molecule), in reasonably good agreement with the measured enthalpies of adsorption [-∆adsH], and with the heterogeneity of surface sites. The possible extrapolation of computed data to vanishing coverage accounts for an initial enthalpy of adsorption larger than 100 kJ/mol. At high coverage (more than one molecule per Ca2+ cation), the BEC of H2O adsorbed on the (001) surface turns out to be larger than for the (010) one and seems to reach a plateau. This was interpreted as due to the onset (at high θ) of lateral interactions among the molecules adsorbed on the (001) surface, which causes the ∆ELC among the adsorbed moieties to contribute significantly to the binding energy (see Table 2 and previous discussion). This does not occur on the (010) surface in which H2O molecules even at high θ follow being engaged with the basic functionalities of the surface. The general agreement between computed BEC and measured [-∆adsH] values is remarkable, in particular when considering that (i) the calorimetric data represent the enthalpy of adsorption of water at 303 K, whereas BEC values are electronic binding energies, which do not include neither zero point energy nor thermal corrections; (ii) B3LYP data have been obtained only on two crystallographic defined surfaces, (001) and (010), whereas the experiments were performed on polycrystalline nanosized samples, for which neither the distribution of different crystal nor the population of possible structural defects are well defined; (iii) the simulated models at high θ are simple microstates on a purely electronic potential energy surface, in that no thermal effects have been included; and (iv) forces deriving from purely London instantaneous dipoles are entirely missing with the B3LYP (and for most of the standard GGA functionals) so that the present BEC are somehow underestimated.

4. Summary and Conclusions The B3LYP computational study of H2O adsorption on hexagonal hydroxyapatite (001) and (010) stoichiometric surfaces has been carried out by means of a periodic approach, using the CRYSTAL06 program. On the (001) surface, the maximum H2O coverage envisaged up to 10 H2O (5 to each slab face) molecules, and, in all cases, H2O remains molecularly adsorbed. The binding energy (in the 65-100 kJ mol-1 per molecule range) derives from the direct interaction of the exposed Ca2+ cations with H2O, which simultaneously engages rather strong H-bonds with the basic oxygen of the PO4 surface group. It is worth noting that at higher

Corno et al.

H2O loadings an important fraction of the binding energy is due to lateral H-bonding interactions among adsorbed H2O molecules. On the (010) surface and for very low coverage (1 H2O per crystal face), H2O chemisorbed on three Ca sites out of the total four available, the Ca1 site being the only surface cation in which H2O is adsorbed molecularly. H2O dissociates without apparent energy barrier and as a highly favored thermodynamic process (the computed reaction energies were in the -320/-250 kJ mol-1 range per H2O molecule adsorbed), yielding new surface terminations, that is, CaOwHw and POHw. The local geometry around the newly formed CaOwHw strongly resembles the environment surrounding the native structural OH group of the (001) surface; that is, the OH is shared between three Ca2+ cations. It is inferred that H2O prefers to stick molecularly on the Ca1 site because this stable structural motif cannot be formed for geometrical reasons. On the reacted (010) surface and at higher θ, H2O is adsorbed molecularly, with binding energy values slightly lower than those computed for the (001) face. Despite the decrease in binding energy at higher H2O loading, the interaction remains favorable on the reacted (010) surface. The reasonable good agreement between the computed and the experimentally measured energetic data on polycrystalline nanosized HA samples allows one to explain the particularly strong affinity of hydroxyapatite for molecular water, as due to the complex bonding mechanism that resulted, at all considered θ, as an interplay between the coordination of H2O onto cus Ca2+ cations through the oxygen end, and a network of H-bonds of variable strength with the basic oxygen atoms of the surface PO4 groups, which cause a reinforcement of the interaction. For both surfaces, IR frequency shifts of the modes of adsorbed H2O with respect to those of free H2O show hypsocromic shifts of the HOH bending mode and large bathocromic shifts of the OH stretching modes due to H-bonds, in semiquantitative agreement with experimental findings. Correlation between OH stretching shifts with the O-H · · · O bond lengths has been shown to hold, as expected from general behavior of H-bond interactions.62 The present study opens up the way for more complex studies of the adsorption of amino acids on the HA surfaces, which is relevant for the interaction between HA and biological systems.63 At the moment, the present results are used in our laboratory to establish the competing adsorption between water and glycine amino acid, with the purpose to assess whether amino acids prefer to interact directly with the HA surface or to be mediated by the molecularly adsorbed water. Acknowledgment. Financial support from the Italian Ministry MIUR (Project COFIN-2006, Prot. 2006032335 Interface phenomena in silica-based nanostructured biocompatible materials contacted with biological systems) and by Regione Piemonte - Italy (Project CIPE-2004, Nanotechnologies and Nanosciences Nanostructured materials biocompatible for biomedical applications) is gratefully acknowledged. Colleagues Gianmario Martra and Luca Bertinetti (Dip. Chimica IFM, University of Torino) are acknowledged for fruitful discussion and for providing the HA samples, synthesized and kindly supplied by ISTEC-CNR (Faenza, Italy). Some of the largest calculations were carried out at the Barcelona Supercomputing Center on the Mare Nostrum resource within project BCV-2008-2-0013: Simulation of peptide folding induced by inorganic materials. LA803253K (63) Rimola, A.; Corno, M.; Zicovich-Wilson, C. M.; Ugliengo, P. J. Am. Chem. Soc. 2008, 130, 16181.