We have carried out “first-principles” Born–Oppenheimer molecular dynamics (BOMD) simulations of the phosphate ions H2PO4– and HPO42– in liquid water and have calculated their IR spectra by Fourier transform techniques from the trajectories. IR bands were assigned by a so-called “generalized normal coordinate analysis”. The effects of including Hartree–Fock (HF) exchange into the density functional theory (DFT) computation of forces were studied by comparing results obtained with the well-known BP, BLYP, and B3LYP functionals. The neglect of dispersion in the functionals was empirically corrected. The inclusion of HF exchange turned out to yield dramatically improved and, thus, quite accurate descriptions of the IR spectra observed for H2PO4– and HPO42– in aqueous solution. An analysis of earlier computational results (Klähn, M. et al. J. Phys. Chem. A 2004, 108, 6186−6194) on these vibrational spectra, which had been obtained in a hybrid setting combining a BP description of the respective phosphate with a simple molecular mechanics (MM) model of its aqueous environment, revealed three different sources of error, (i) the BP force field of the phosphates is much too soft and would have required a substantial scaling of frequencies, (ii) the oversimplified water force field entailed incorrect solvation structures and, thus, qualitatively wrong patterns of solvatochromic band shifts, and (iii) quantitative frequency computations additionally required the inclusion of HF exchange. Thus, the results of the B3LYP BOMD simulations do not only characterize physical properties like the IR spectra or the solvation structures of the phosphate systems but also provide clues for the future design of simplified but nevertheless reasonably accurate DFT/MM methods applicable to phosphates.