Periodic plane-wave density functional theory (DFT) and molecular cluster hybrid molecular orbital-DFT (MO-DFT) calculations were performed on models of phosphate surface complexes on the (100), (010), (001), (101), and (210) surfaces of alpha-FeOOH (goethite). Binding energies of monodentate and bidentate HPO42- surface complexes were compared to H2PO4- outer-sphere complexes. Both the average potential energies from DFT molecular dynamics (DFT-MD) simulations and energy minimizations were used to estimate adsorption energies for each configuration. Molecular clusters were extracted from the energy-minimized structures of the periodic systems and subjected to energy reminimization and frequency analysis with MO-DFT. The modeled P-O and P---Fe distances were consistent with EXAFS data for the arsenate oxyanion that is an analog of phosphate, and the interatomic distances predicted by the clusters were similar to those of the periodic models. Calculated vibrational frequencies from these clusters were then correlated with observed infrared bands. Configurations that resulted in favorable adsorption energies were also found to produce theoretical vibrational frequencies that correlated well with experiment. The relative stability of monodentate versus bidentate configurations was a function of the goethite surface under consideration. Overall, our results show that phosphate adsorption onto goethite occurs as a variety of surface complexes depending on the habit of the mineral (i.e., surfaces present) and solution pH. Previous IR spectroscopic studies may have been difficult to interpret because the observed spectra averaged the structural properties of three or more configurations on any given sample with multiple surfaces.