  Octahedral Fe$^{2+}$ molecules are particularly interesting as they often
exhibit a spin-crossover transition. In spite of the many efforts aimed at
assessing the performances of density functional theory for such systems, an
exchange-correlation functional able to account accurately for the energetic of
the various possible spin-states has not been identified yet. Here we
critically discuss the issues related to the theoretical description of this
class of molecules from first principles. In particular we present a comparison
between different density functionals for four ions, namely
[Fe(H$_2$O)$_6$]$^{2+}$, [Fe(NH$_3$)$_6$]$^{2+}$, [Fe(NCH)$_6$]$^{2+}$ and
[Fe(CO)$_6$]$^{2+}$. These are characterized by different ligand-field
splittings and ground state spin multiplicities. Since no experimental data are
available for the gas phase, the density functional theory results are
benchmarked against those obtained with diffusion Monte Carlo, one of the most
accurate methods available to compute ground state total energies of quantum
systems. On the one hand, we show that most of the functionals considered
provide a good description of the geometry and of the shape of the potential
energy surfaces. On the other hand, the same functionals fail badly in
predicting the energy differences between the various spin states. In the case
of [Fe(H$_2$O)$_6$]$^{2+}$, [Fe(NH$_3$)$_6$]$^{2+}$, [Fe(NCH)$_6$]$^{2+}$, this
failure is related to the drastic underestimation of the exchange energy.
Therefore quite accurate results can be achieved with hybrid functionals
including about 50% of Hartree-Fock exchange. In contrast, in the case of
[Fe(CO)$_6$]$^{2+}$, the failure is likely to be caused by the
multiconfigurational character of the ground state wave-function and no
suitable exchange and correlation functional has been identified.
