A method based on experimental and in silico evaluations for investigating interactions of organic phosphates and phosphonates with hydroxyapatite was developed. This quick and easy method is used for determination of differences among organophosphorus compounds of various structures in their mineral binding affinities. Empirical sorption evaluation was carried out using liquid chromatography with tandem mass spectrometry or UV–VIS spectroscopy. Raman spectroscopy was used to confirm sorption of organic phosphates and phosphonates on hydroxyapatite. Polymer-ceramic monolithic material and bulk hydroxyapatite were applied as sorbent materials. Furthermore, a Polymer-ceramic Monolithic In-Needle Extraction device was used to investigate both sorption and desorption steps. Binding energies were computed from the fully optimised structures utilising Density Functional Theory (DFT) at B3LYP/6-31 + G(d,p) level. Potential pharmacologic and toxic effects of the tested compounds were estimated by the Prediction of the Activity Spectra of Substances using GeneXplain software.