The influence of ligands on electronic structure of small gold clusters (Au₂, Au₄) has been investigated by density functional theory (DFT). Specifically, we study the effect of bonding of four donor ligands (NH₃, NMe₃, PH₃, and PMe₃) on cluster geometries and energetics in gas phase and in solution. Performance of five generations of DFT functionals and five different basis sets is assessed. Our results benchmark the importance of the DFT functional model and polarization functions in the basis set for calculations of ligated gold cluster systems. We obtain NMe₃ ≈ NH₃ < PH₃ < PMe₃ order of ligand binding energies and observe shallow potential energy surfaces in all molecules. The latter is likely to lead to a conformational freedom in larger clusters with many ligands in solution at ambient conditions. The study suggests appropriate quantum-chemical methodology to reliably model small noble metal clusters in a realistic ligand environment typically present in experiments.