We present a non-orthogonal fragment ab initio methodology for the calculation of crystal field energy levels and magnetic properties in lanthanide complexes, implementing a systematic description of non-covalent contributions to metal-ligand bonding. The approach has two steps. In the first step, appropriate ab initio wavefunctions for the various ionic fragments (lanthanide ion and coordinating ligands) are separately optimized, accounting for the electrostatic influence of the surrounding environment, within various approximations. In the second and final step, the scalar relativistic (DKH2) electrostatic Hamiltonian of the whole molecule is represented on the basis of the optimized metal-ligand multiconfigurational non-orthogonal group functions (MC-NOGF), and reduced to an effective (2J+1)-dimensional non-orthogonal Configuration Interaction (CI) problem via L{\"o}wdin-partitioning. Within the proposed formalism, the projected Hamiltonian can be implemented to any desired order of perturbation theory in the fragment-localised excitations out of the degenerate space, and its eigenvalues and eigenfunctions are systematic approximations to the crystal field energies and wavefunctions. We present a preliminary implementation of the proposed MC-NOGF method to first-order degenerate perturbation theory within our own ab initio code CERES, and compare its performance both with the simpler non-covalent orthogonal ab initio approach Fragment Ab Initio Model Potential (FAIMP) approximation, and with the full CAHF/CASCI-SO method, accounting for metal-ligand covalency in a mean-field manner. We find that energies and magnetic properties for 44 complexes obtained via an iteratively optimized version of our MC-NOGF first-order non-covalent method, compare remarkably well to the full CAHF/CASCI-SO method including metal-ligand covalency, and are superior to the best purely electrostatic results achieved via an iteratively optimized version of the FAIMP approach.<br>