We have performed density functional theory (DFT) calculations of iron-porphyrin (FeP) and its complexes with O-2, CO, NO, and imidazole(Im). Our fully optimized structures agree well with the available experimental data for synthetic heme models. Comparison with crystallographic data for proteins highlights interesting features of carbon monoxymyoglobin. The diatomic molecule induces a 0.3-0.4 Angstrom displacement of the Fe atom out of the porphyrin nitrogen (N-p) plane and a doming of the overall porphyrin ring. The energy of the iron-diatomic bond increases in the order Fe-O-2 (9 kcal/mol) < Fe-CO (26 kcal/mol) < Fe-NO (35 kcal/mol). The ground state of FeP(O-2) is an open shell singlet. The bent Fe-O-2 bond can be formally described as Fe-III-O-2(-), and it is characterized by the anti-ferromagnetic coupling between one of the d electrons of Fe and one of the pi* electrons of O-2 FeP(CO) is a closed shell singlet, with a linear Fe-C-Q bond, The complex with NO has a doublet ground state and a Fe-NO geometry intermediate between that of FeP(CO) and FeP(O-2). The bending of the Fe-(diatomic) angle requires a rather low energy for these three complexes, resulting in large-amplitude oscillations of the ligand even at room temperature. The addition of an imidazole ligand to FeP moves the Fe atom out of the porphyrin plane toward the imidazole and decreases significantly the energy differences among the spin states. Moreover, our calculations underline the potential role of the imidazole ligand in controlling the electronic structure of FeP by changing the out-of-planarity of the Fe atom. The presence of the imidazole increases the strength of the Fe-O-2 and Fe-CO bonds (15 and 35 kcal/mol, respectively), but does not affect the energy of the Fe-NO bond, while the resulting FeP(Im)(NO) complex exhibits a longer and weaker Fe-Im bond.