Understanding the structures of phosphate glasses is important to many of their technological applications. Molecular dynamics simulations are commonly used to generate structure models of sodium phosphate glasses, and those with partial charge pairwise potentials have been successfully applied for modeling other network glasses, such as silicate and aluminosilicate glasses. In this work, we show that the addition of a three-body term is essential in regulating the intertetrahedral bond angles, as well as Qn speciation in comparison to experiments. Simulation results with and without three-body terms were compared and validated with experimental results, including neutron structure factors. Further comparison with glass structures fully relaxed with first-principles density functional theory was performed to evaluate the simulation results. The results show that the addition of three-body terms is vital for the modeling of phosphate glasses, and it can significantly improve the description of short- and medium-range structures and properties.