Recently, an alternative strategy has been proposed by us , 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 . The optimizations of the MO and CI coefficients are then decoupled, and the latter coefficients are determined in a second step by direct diagonalization of the total Hamiltonian within the Complete Active Space Configuration In-teraction with Spin-Orbit coupling (CASCI-SO) method. CAHF/CASCI-SO method has been shown to produce consistent results with CASSCF/RASSI-SO approach with a significant saving in computational time, and it has been implemented into the Computational Emulator of Rare Earth Systems (CERES) quantum chemistry package , an ad hoc tool specifically designed for the calculation of the electronic structure and magnetic properties of lanthanide complexes.