Molecular dynamics simulations of molten hydroxyapatite were performed, for the first time, in the range 2000 K < T < 3000 K and pres-sures up to 20 GPa. The all-atom Born–Huggins–Mayer potential energy function employed had been previously used to study the ther-modynamic 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, α p and the isothermal compressibility, κ T , could be evaluated. The heat capacity at room pressure, C p , in the range 2000–3000 K, was estimated from the plot of the molar enthalpy of the melt as a function of temperature, H m = A 0 + AT + BT 2 + C/T (A 0 = −3.7490 × 10 4 kJ mol −1 , A = 3.5842 kJ mol −1 K −1 , B = −5.6989 × 10 −4 kJ mol −1 K −2 , C = −3.0061 × 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 ρ 2 , has been established for molten hydroxyapatite.