Phosphate chemistry is involved in many key biological processes yet the underlying mechanism often remains unclear. For theoretical analysis to effectively complement experimental mechanistic analysis, it is essential to develop computational methods that can capture the complexity of the underlying potential energy surface and allow for sufficient sampling of the configurational space. To this end, we report the parameterization of an approximate density functional theory, Self-Consistent-Charge Density-Functional Tight-Binding (SCC-DFTB) method for systems containing phosphorus. Compared to high-level density functional theory and ab initio (MP2 and G3B3) results, the standard second-order parameterization is shown to give reliable structures for a diverse set of phosphate compounds but inaccurate energetics. With the on-site third-order terms included, referred to as SCC-DFTBPA, calculated proton affinities of phosphate compounds are substantially improved, although it remains difficult to obtain reliable proton affinity for both phosphates and compounds that do not contain phosphorus, indicating that further improvement in the formulation of SCC-DFTB is still a challenge to meet. To make SCC-DFTB applicable to phosphate reactions in the current (on-site-third-order-only) formulation, a "reaction-specific" parameterization, referred to as SCC-DFTBPR, is developed based on hydrolysis reactions of model phosphate species. Benchmark calculations in both the gas-phase and solution-phase indicate that SCC-DFTBPR gives reliable structural properties and semi-quantitative energetics for phosphate hydrolysis reactions. Since the number of reaction-specific parameters is small, it is likely that SCC-DFTBPR is applicable to a broad set of phosphate species. Indeed, for 56 reaction exothermicities and 47 energy barriers related to RNA catalysis model reactions collected from the QCRNA database, which involve molecules rather different from those used to parameterize SCC-DFTBPR, the corresponding root-mean-square difference between SCC-DFTBPR and high-level DFT results is only 5.3 kcal/mol. We hope that the parameterized SCC-DFTB models will complement NDDO based reaction-specific models (e.g., AM1-d/PhoT) and high-level ab initio QM/MM methods in better understanding the mechanism of phosphate chemistry in condensed phase, particularly biological systems.