The performance of density functional theory (DFT) has been widely examined with regard to its ability to predict the properties of minerals, though less attention has been given to the correct determination of the relative stability of structurally similar polymorphs. Here a detailed examination is performed of the numerical and theoretical factors that may influence the structure and relative energetics of two such polymorphs of iron disulfide, namely pyrite and marcasite, within density functional theory. Both the local density approximation and commonly used generalized gradient approximation exchange-correlation functionals, such as PBE, are found to predict that marcasite is more stable than pyrite, at variance with experiment. Allowing for the zero point energy of vibration fails to remedy this discrepancy. While inclusion of a 2 sufficiently large Hubbard U parameter for iron is found to reverse the stability, this comes at the expense of a very poor description of other properties. Examination of three generalized gradient approximations developed specifically for the solid-state, namely AM05, Wu-Cohen and PBEsol, demonstrates that all of these functionals offer a superior description of the structures and relative energies of pyrite and marcasite through correctly predicting that the former is the ground state phase at ambient conditions.