Molecular dynamics simulations of molten hydroxyapatite were performed, for the first time, in the range 2000 K < T < 3000 K and pressures up to 20 GPa. The all-atom Born-Huggins-Mayer potential energy function employed had been previously used to study the thermodynamic properties of the solid compound. High-temperature simulation runs were used to generate the p-V-m-T surface of the melt, from which properties like the isobaric thermal expansion coefficient, alpha(p) and the isothermal compressibility, K-T, could be evaluated. The heat capacity at room pressure, C-p, in the range 2000-3000K, was estimated from the plot of the molar enthalpy of the melt as a function of temperature, H-m =A(0) + AT + BT2 + C/T (A(0) = -3.7490 x 10(4) kJ mol(-1). A = 3.5842 kJ mol(-1) K-1. B = -5.6989 x 10(-4) kJ mol(-1) K-2, C = -3.0061 x 10(5) kJ mol(-1) K). C-p varies from 1373 J mol(-1) K-1 (T = 2000 K) to 180 J mol(-1) K-1 (T = 3000 K). The intermolecular atom-atom distribution functions, at several temperatures and pressures, were also investigated. A universal EoS proposed by Parsafar et al. was shown to give a good account of the MD data, the precision being better than 0.5%. Likewise, the Parsafar-Mason regularity which assumes a linear dependence of (Z-1) V-2 on rho(2), has been established for molten hydroxyapatite. (c) 2005 Elsevier B.V. All rights reserved.