A theoretical, computational, and conceptual framework for the interpretation and prediction of the magnetic anisotropy of transition metal complexes with orbitally degenerate or orbitally nearly degenerate ground states is explored. The treatment is based on complete active space self-consistent field (CASSCF) wave functions in conjunction with N-electron valence perturbation theory (NEVPT2) and quasidegenerate perturbation theory (QDPT) for treatment of magnetic field- and spin-dependent relativistic effects. The methodology is applied to a series of Fe(II) complexes in ligand fields of almost trigonal pyramidal symmetry as provided by several variants of the tris-pyrrolylmethyl amine ligand (tpa). These systems have recently attracted much attention as mononuclear single-molecule magnet (SMM) complexes. This study aims to establish how the ligand field can be fine tuned in order to maximize the magnetic anisotropy barrier. In trigonal ligand fields high-spin Fe(II) complexes adopt an orbitally degenerate E-5 ground state with strong in-state spin-orbit coupling (SOC). We study the competing effects of SOC and the E-5 circle times epsilon multimode Jahn-Teller effect as a function of the peripheral substituents on the tpa ligand. These subtle distortions were found to have a significant effect on the magnetic anisotropy. Using a rigorous treatment of all spin multiplets arising from the triplet and quintet states in the d(6) configuration the parameters of the effective spin-Hamiltonian (SH) approach were predicted from first principles. Being based on a nonperturbative approach we investigate under which conditions the SH approach is valid and what terms need to be retained. It is demonstrated that already tiny geometric distortions observed in the crystal structures of four structurally and magnetically well-documented systems, reported recently, i.e., [Fe(tpa(R))](-) (R = tert-butyl, Tbu (1), mesityl, Mes (2), phenyl, Ph (3), and 2,6-difluorophenyl, Dfp (4), are enough to lead to five lowest and thermally accessible spin sublevels described sufficiently well by S = 2 SH provided that it is extended with one fourth order anisotropy term. Using this most elementary parametrization that is consistent with the actual physics, the reported magnetization data for the target systems were reinterpreted and found to be in good agreement with the ab initio results. The multiplet energies from the ab initio calculations have been fitted with remarkable consistency using a ligand field (angular overlap) model (ab initio ligand field, AILFT). This allows for determination of bonding parameters and quantitatively demonstrates the correlation between increasingly negative D values and changes in the sigma-bond strength induced by the peripheral ligands. In fact, the sigma-bonding capacity (and hence the Lewis basicity) of the ligand decreases along the series 1 > 2 > 3 > 4.