Due to the typical over-parametrization afflicting phenomenological crystal field theory models , nowadays the ab initio cal-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 . 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 , and it has extensively used to simulate the magnetic properties of lanthanide-based systems.