This paper reports a new AM1/d model for phosphorus that can be used to model nucleophilic attack of phosphates relevant for biological phosphate hydrolysis reactions. The parameters were derived from a quantum dataset calculated with hybrid density-functional theory [B3LYP/6-311++G(3df,2p)//B3LYP/6-31++G(d,p)] of phosphates and phosphoranes in various charge states, and on transitions states for nucleophilic attacks. A suite of non-linear optimization methods is outlined for semiempirical parameter development based on integrated evolutionary (genetic), Monte Carlo simulated annealing and direction set minimization algorithms. The performance of the new AM1/d model and the standard AM1 and MNDO/d models are compared with the density-functional results. The results demonstrate that the strategy of developing semiempirical parameters specific for biological reactions offers considerable promise for application to large-scale biological problems.