Equations are derived and discussed that allow the computation of zero-field splitting (ZFS) tensors in transition metal complexes for any value of the ground-state total spin S. An effective Hamiltonian technique is used and the calculation is carried to second order for orbitally nondegenerate ground states. The theory includes contributions from excited states of spin S and S +/- 1. This makes the theory more general than earlier treatments. Explicit equations are derived for the case where all states are well described by single-determinantal wave functions, for example restricted open shell Hartree-Fock (HF) and spin-polarized HF or density functional (DFT) calculation schemes. Matrix elements are evaluated for many electron wave functions that result from a molecular orbital (MO) treatment including configuration interaction (CI). A computational implementation in terms of bonded functions is outlined. The problem of ZFS in high-spin ferric complexes is treated at some length, and contributions due to low-symmetry distortions, anisotropic covalency, charge-transfer states, and ligand spin-orbit coupling are discussed. ROHF-INDO/S-Cl results are presented for FeCl4- and used to evaluate the importance of the various terms. Finally, contributions to the experimentally observed reduction of the metal spin-orbit coupling constants (the relativistic nephelauxetic effect) are discussed. B3LYP and Hartree-Fock calculations for FeCl4- are used to characterize the change of the iron 3d radial function upon complex formation. It is found that the iron 3d radial distribution function is significantly expanded and that the expansion is anisotropic. This is interpreted as a combination of reduction in effective charge on the metal 3d electrons (central field covalence) together with expansive promotion effects that are a necessary consequence of chemical bond formation. The [r(-3)](3d) values that are important in the interpretation of magnetic data are up to 15% reduced from their free-ion value before any metal-ligand orbital mixing (symmetry-restricted covalency) is taken into account. Thus the use of free-ion values for spin-orbit coupling and related constants in the analysis of experimental data leads to values for MO coefficients that overestimate the metal-ligand covalency.