We present a number of computationally cost-effective approaches to calculate magnetic excitations (i.e. crystal field energies and magnetic anisotropies in the lowest spin-orbit multiplet) in lanthanide complexes. In particular, we focus on the representation of the spin-orbit coupling term of the molecular Hamiltonian, which has been implemented within the quantum chemistry package CERES using various approximations to the Breit-Pauli Hamiltonian. The approximations include the (i) bare one-electron approximation, (ii) atomic mean field and molecular mean field approximations of the two-electron term, (iii) full representation of the Breit-Pauli Hamiltonian. Within the framework of the CERES implementation, the spin-orbit Hamiltonian is always fully diagonalized together with the electron repulsion Hamiltonian (CASCI-SO) on the full basis of Slater determinants arising within the 4 f ligand field space. For the first time, we make full use of the Cholesky decomposition of two-electron spin-orbit integrals to speed up the calculation of the two-electron spin-orbit operator. We perform an extensive comparison of the different approximations on a set of lanthanide complexes varying both the lanthanide ion and the ligands. Surprisingly, while our results confirm the need of at least a mean field approach to accurately describe the spin-orbit coupling interaction within the ground Russell-Saunders term, we find that the simple bare one-electron spin-orbit Hamiltonian performs reasonably well to describe the crystal field split energies and g tensors within the ground spin-orbit multiplet, which characterize all the magnetic excitations responsible for lanthanide-based single-molecule magnetism. Electronic Supplementary Information (ESI) available: [details of any supplementary information available should be included here]. See DOI: 00.0000/00000000. culations of the crystal field energy levels of lanthanide complexes play an important role in the interpretation of the experimental results. Among the various available ab initio approaches, the Complete Active Space Self-Consistent Field with Restricted Active Space State Interaction via Spin-Orbit coupling (CASSCF/RASSI-SO) method has been successfully used in the literature to compute the magnetic properties of the lanthanide complexes 9-11 . Such a method involves converging a set of molecular orbitals (MO) for each spin configuration in the CASSCF step. Then diagonalizing the Spin-Orbit Coupling (SOC) operator over the basis of the CASSCF wavefunctions. This is the method of choice in the popular package MOLCAS 12 , and it has extensively used to simulate the magnetic properties of lanthanide-based systems.Recently, an alternative strategy has been proposed by us 13,14 , were a set of MO are optimized by minimizing the average-energy functional represented on the basis of all the possible CI configurations, which is called Configuration Averaged Hartree Fock (CAHF) method 13,15,16 . The optimizations of the MO and CI coefficients are then decoupled, and the latte...