Due to the importance of the understanding of dissolution behavior of phosphate-based bioglasses (PBGs) in different biomedical applications, binary sodium and calcium phosphate glasses have been simulated for the first time using a newly developed ReaxFF force field and a standard melt-quench method with the LAMMPS classical molecular dynamics software. The partial radial distribution function of P-O within the first coordination shell indicated two distinct peaks corresponding to phosphorus bonding to NBO and BO, respectively, at distances consistent with those observed experimentally and a P-O coordination number of 4.0. Angular distribution functions were consistent with the experimental data. The calculated network connectivities are in good agreement with experimental data, and the detailed Qn distributions are broader than would be expected.