Antiferromagnetic coupling in multinuclear transition metal complexes usually leads to electronic ground states that cannot be described by a single Slater determinant and that are therefore difficult to describe by Kohn-Sham density functional methods. Density functional calculations in such cases are usually converged to broken symmetry solutions which break spin and, in many cases, also spatial symmetry. While a procedure exists to extract isotropic Heisenberg (exchange) coupling constants from such calculations, no such approach is yet established for the calculation of magnetic anisotropy energies or zero field splitting parameters. This work proposes such a procedure. The broken symmetry solutions are not only used to extract the exchange couplings but also single-ion D tensors which are then used to construct a (phenomenological) spin Hamiltonian, from which the magnetic anisotropy and the zero-field energy levels can be computed. The procedure is demonstrated for a bi- and a trinuclear Mn(III) model compound.