A Molecular Dynamics Study of the Thermodynamic Properties of Calcium Apatites. 1. Hexagonal Phases Fernando J. A. L. Cruz, Jose ´ N. Canongia Lopes, and Jorge C. G. Calado* Centro de Quı ´mica Estrutural, Complexo Interdisciplinar, Instituto Superior Te ´ cnico, 1049-001 Lisboa, Portugal Manuel E. Minas da Piedade Departamento de Quı ´mica e Bioquı ´mica, Faculdade de Cie ˆ ncias, UniVersidade de Lisboa, 1649-016 Lisboa, Portugal ReceiVed: August 3, 2005; In Final Form: September 30, 2005 Structural and thermodynamic properties of crystal hexagonal calcium apatites, Ca 10 (PO 4 ) 6 (X) 2 (X ) OH, F, Cl, Br), were investigated using an all-atom Born-Huggins-Mayer potential by a molecular dynamics technique. The accuracy of the model at room temperature and atmospheric pressure was checked against crystal structural data, with maximum deviations of ca. 4% for the haloapatites and 8% for hydroxyapatite. The standard molar lattice enthalpy, Δ lat H 298 °, of the apatites was calculated and compared with previously published experimental results, the agreement being better than 2%. The molar heat capacity at constant pressure, C p,m , in the range 298-1298 K, was estimated from the plot of the molar enthalpy of the crystal as a function of temperature, H m ) (H m,298 - 298C p,m ) + C p,m T, yielding C p,m ) 694 ( 68 J‚mol -1 ‚K -1 , C p,m ) 646 ( 26 J‚mol -1 ‚K -1 , C p,m ) 530 ( 34 J‚mol -1 ‚K -1 , and C p,m ) 811 ( 42 J‚mol -1 ‚K -1 for hydroxy-, fluor-, chlor-, and bromapatite, respectively. High-pressure simulation runs, in the range 0.5-75 kbar, were performed in order to estimate the isothermal compressibility coefficient, κ T , of those compounds. The deformation of the compressed solids is always elastically anisotropic, with BrAp exhibiting a markedly different behavior from those displayed by HOAp and ClAp. High-pressure p-V data were fitted to the Parsafar-Mason equation of state with an accuracy better than 1%. 1. Introduction Apatites are a class of compounds of general formula M 10 (PO 4 ) 6 X 2 , where M is a divalent metal cation (primarily Ca 2+ , Sr 2+ , Ba 2+ , Cd 2+ , or Pb 2+ ), and X is a monovalent anion (OH - ,F - , Cl - , and Br - ). The calcium derivatives are by far the most abundant, and they constitute one of the principal host phases for condensed phosphorus and halogens in the solar system, 1-4 at temperatures below ca. 700 K. Hydroxyapatite (HOAp) is a major constituent of mammalian bone and teeth enamel, where it occurs as a nonstoichiometric, amorphous, and carbonated form, but it is also very important as a biomaterial for medical implants. 5-8 Fluorapatite (FAp) is the most com- monly occurring calcium apatite, whereas both pure chlorapatite (ClAp) and bromapatite (BrAp) do not occur naturally on Earth, although they can be found in some astral bodies such as meteorites. 3 The structure of the apatites is basically a crystal lattice composed of tetrahedral phosphate anions and calcium cations, with the hydroxyl and halide ions lying on the crystallographic c-axis perpendicular to the planes defined by three calcium ions. The stacking of these planes results in a series of parallel hexagonal pipes where the substituents are located (Figure 1a). While the fluoride ions are exactly located on the triangle centers, the hydroxyl, chloride, and bromide ions are increasingly shifted toward a position above those triangles (Figure 1b), resulting in an overall hexagonal symmetry (P6 3 / m). As far as the authors are aware, only hydroxy- and fluorapatite have been the subject of molecular dynamics studies. On the basis of structural data and electronic structure calculations, Mkhonto 9 and de Leeuw 10 have proposed potential energy models for FAp and HOAp, respectively. Hauptmann et al. 11 studied both monoclinic hydroxyapatite and hexagonal fluora- patite and fitted experimental data and quantum-chemical calculation results to describe the solids using a Born- Huggins-Mayer type potential. Meis et al. 12 developed a similar potential energy function for fluorapatite, modeling the phos- phate tetrahedra with a three-body short-range harmonic po- tential. The purpose of this work is to calculate the cohesive energy of the apatites and study the evolution of their structures under high-temperature or high-pressure conditions by assuming the Born model of ionic solids 13 and using a molecular dynamics (MD) technique to perform the relevant calculations. Simulation runs at room temperature (T ) 298 K) and atmospheric pressure (p ) 1 bar) were used to calculate the standard molar lattice enthalpy of the different crystals; high-temperature simulations in the range 298 K < T < 1298 K were used to estimate their heat capacity at constant pressure; and high-pressure runs in the range 0.5 kbar < p < 75 kbar were performed to yield the corresponding isothermal compressibility values. 2. Theoretical Methods Force Field Model. All apatite molecules in this work were modeled by an all-atom force-field based on a four-parameter * To whom correspondence may be addressed. E-mail: jcalado@ist.utl.pt. 24473 J. Phys. Chem. B 2005, 109, 24473-24479 10.1021/jp054304p CCC: $30.25 © 2005 American Chemical Society Published on Web 12/02/2005