Previous work testing density functionals for use in calculating high-spin-low-spin energy differences, deltaE(HL), for iron(II) spin-crossover transitions has tended to conclude that only properly reparametrized hybrid functionals can predict deltaE(HL) since it seems to depend critically on a correct description of the electron pairing energy governed by the exchange term. Exceptions to this rule are the previous three papers (I, II, and III in the present series of papers) where it was found that modern generalized gradient approximations (GGAs) and meta-GGAs could do as well as hybrid functionals, if not better, for this type of problem. In the present paper, we extend these previous studies to five more molecules which are too large to treat with high-quality ab initio calculations, namely, the series [Fe(L)('NHS(4)')], where NHS(4)=2.2'-bis(2-mercaptophenylthio)diethylamine dianion, and L=NH(3), N(2)H(4), PMe(3), CO, and NO(+). Since we know of no reliable experimental estimate of deltaE(HL), we content ourselves with a comparison against the experimentally determined ground-state spin symmetry including, in so far as possible, finite-temperature effects. Together with the results of Papers I, II, and III, this paper provides a test of a large number of functionals against the high-spin/low-spin properties of a diverse set of Fe(II) compounds, making it possible to draw some particularly interesting conclusions. Trends among different classes of functionals are discussed and it is pointed out that there is at least one functional, namely, the OLYP generalized gradient approximation, which is able to give a reasonably good description of the delicate spin energetics of Fe(II) coordination compounds without resorting to hybrid functionals which require the relatively more expensive calculation of a Hartree-Fock-type exchange term.