Selected theoretical methods, basis sets and solvation models have been tested in their ability to predict 31P NMR chemical shifts of large phosphorous-containing molecular systems in solution. The most efficient strategy was found to involve NMR shift calculations at the GIAO-MPW1K/6-311++G(2d,2p)//MPW1K/6-31G(d) level in combination with a dual solvation model including the explicit consideration of single solvent molecules and a continuum (PCM) solvation model. For larger systems it has also been established that reliable 31P shift predictions require Boltzmann averaging over all accessible conformations in solution.