We report modeling studies of high spin and low spin forms of the test set of problematic cases of pseudo-octahedral Fe(II) complexes defined by Ye and Neese. These include the most thoroughly studied spin crossover compound Fe(SCN)(2)(Phen)(2) and a more recently discovered Fe(II)-2,13-dimethyl-6,9-dio xa-3,12,18-triaza-bicyclo[12.3.1]octadeca-1(18),2,12,14,16-pentene. Reliable estimates of the often small high-spin to low-spin energy gap require balanced treatment of electronic structure in very different bonding environments. Density functionals differ in their bias toward low spin or high spin configurations. For the test set defined by Ye and Neese, we describe the estimates of their energy gaps between high-spin and low-spin forms for 75 functionals; geometries are obtained by optimization with B3LYP with zero-order relativistic approximation (ZORA) in a triple-zeta quality basis. Fairly reliable results for high-spin to low-spin energy gaps are obtained by functionals IMP, HCTH, OPBE, and B3LYP*, and with the omega B97XD functional expressed in the cc-pVTZ basis, about half the size of the def2-QZVPP used by Ye and Neese. The best functionals are still biased systematically toward the high spin by about 20 kJ/mol. We apply these methods in a new study of a pseudo C-3v four coordinate crossover species reported by Scepaniak et al. Their four coordinate spin crossover complex of Fe with aryl-substituted scorpionate and N = P(Ph)(3) ligands is modeled by a truncated structure lacking aromatic substituents. These substituents appear to stabilize the lower spin. Reproducing the 80 K crossover temperature requires an empirical correction of about 12 kJ/mot in favor of the singlet for this truncated model system. (C) 2015 Elsevier B.V. All rights reserved.