Electron spins in molecules and materials may be manipulated and used to store information, and hence they are interesting resources for quantum technologies. A way to understand the physical properties of electron spins is to probe their interactions with electromagnetic fields. Such interactions can be described by using a so-called spin Hamiltonian, with parameters derived from either experiments or calculations. For a single electron spin (e.g., associated to a point-defect in a semiconductor or insulator), the leading terms in the spin Hamiltonian arewhere µ B is the Bohr magneton, S is the electron spin operator, B is an external magnetic field, g and D are rank-2 tensors that characterize the strength of the Zeeman interaction, and the zero-field splitting (ZFS), respectively. Experimentally, the spin Hamiltonian parameters g and D may be obtained by electron paramagnetic resonance (EPR). The ZFS tensor describes the lifting of degeneracy of spin sublevels in the absence of external magnetic fields, and is an important property of open-shell molecules and spin defects in semiconductors with spin quantum number S ≥ 1. The ZFS tensor can be predicted from first-principles calculations, thus complementing experiments and providing valuable insight into the design of novel molecules and materials with desired spin properties. Furthermore, the comparison of computed and measured ZFS tensors may provide important information on the atomistic structure and charge state of defects in solids, thus helping to identify the defect configuration present in experimental samples. Therefore, the development of robust methods for the calculation of the ZFS tensor is an interesting topic in molecular chemistry and materials science.In this work we describe the code PyZFS for the calculation of the ZFS tensor D of molecules and solids, based on wavefunctions obtained from density functional theory (DFT) calculations. For systems without heavy elements, i.e., where spin-orbit coupling is negligible, magnetic spinspin interactions are the dominant ones in the determination of the ZFS tensor. For molecules and materials with magnetic permeability close to the vacuum permeability µ 0 , the spin-spin ZFS tensor evaluated using the DFT Kohn-Sham wavefunctions, is given by Harriman (1978):where a, b = x, y, z are Cartesian indices; γ e is the gyromagnetic ratio of electrons; the summation runs over all pairs of occupied Kohn-Sham orbitals; χ ij = ±1 for parallel and antiparallel spins, respectively; Φ ij (r, r ′ ) are 2 × 2 determinants formed from Kohn-Sham orbitals ϕ i and ϕ j , Φ ij (r, r ′ ) = 1 √ 2 [ ϕ i (r)ϕ j (r ′ ) − ϕ i (r ′ )ϕ j (r)].Ma et al., (2020). PyZFS: A Python package for first-principles calculations of zero-field splitting tensors.