We report the results of a series of density functional theory (DFT) calculations aimed at predicting the Fe-57 Mossbauer electric field gradient (EFG) tensors (quadrupole splittings and asymmetry parameters) and their orientations in S = 0, 1/2, 1, 3/2, 2, and 5/2 metalloproteins and/or model systems. Excellent results were found by using a Wachter's all electron basis set for iron, 6-311G* for other heavy atoms, and 6-31G* for hydrogen atoms, BPW91 and B3LYP exchange-correlation functionals, and spin-unrestricted methods for the paramagnetic systems. For the theory versus experiment correlation, we found R-2 = 0.975, slope = 0.99, intercept = -0.08 mm sec(-1), rmsd = 0.30 mm sec(-1) (N = 23 points) covering a DeltaE(Q) range of 5.63 mm s(-1) when using the BPW91 functional and R-2 = 0.978, slope = 1.12, intercept = -0.26 mm sec(-1), rmsd = 0.31 mm sec(-1) when using the B3LYP functional. DeltaE(Q) values in the following systems were successfully predicted: (1) ferric low-spin (S = 1/2) systems, including one iron porphyrin with the usual (d(xy))(2)(d(xy)d(yz))(3) electronic configuration and two iron porphyrins with the more unusual (d(xz)d(yz))(4)(d(xy))(1) electronic configuration; (2) ferrous NO-heme model compounds (S = 1/2); (3) ferrous intermediate spin (S = 1) tetraphenylporphinato iron(II); (4) a ferric intermediate spin (S = 3/2) iron porphyrin; (5) ferrous high-spin (S = 2) deoxymyoglobin and deoxyhemoglobin; and (6) ferric high spin (S = 5/2) metmyoglobin plus two five-coordinate and one six-coordinate iron porphyrins. In addition, seven diamagnetic (S = 0, d(6) and d(8)) systems studied previously were reinvestigated using the same functionals and basis set scheme as used for the paramagnetic systems. All computed asymmetry parameters were found to be in good agreement with the available experimental data as were the electric field gradient tensor orientations. In addition, we investigated the electronic structures of several systems, including the (dxy)(2) (d(xz),d(yz))(3) and (d(xz),d(yz))(4)(d(xy))(1) [Fe(III)/porphyrinate](+) cations as well as the NO adduct of Fe(II)(octaethylporphinate), where interesting information on the spin density distributions can be readily obtained from the computed wave functions.