The energy levels of Ln³⁺ ions are known to be only slightly dependent on the ion environment. This allows one to predict the spectra of f–f transitions in Ln³⁺ complexes using group theory and simple semiempirical models: Russell–Saunders scheme for spin–orbit coupling, ligand-field theory for the splitting of the electronic levels, and Judd–Ofelt parameterization for reproducing the intensity of f–f transitions. Nevertheless, a fully ab initio computational scheme employing no empirical parameterization and suitable for any asymmetrical environment of Ln³⁺ would be instructive. Here we present such a scheme based on the multireference SA-CASSCF/XMCQPDT2/SO-CASSCF (state-averaged complete active space SCF, quasi-degenerate perturbation theory, and spin–orbit CASSCF) approach for trivalent lanthanide ions from Ce³⁺ (4f¹) to Yb³⁺ (4f¹³). To achieve the most accurate results, we analyse the factors that influence the accuracy of the calculation: basis set size, state averaging scheme, effect of the low-spin states on the energy gap between the high-spin states (e.g., effect of triplets on the septet–quintet gaps in f⁶ or f⁸ configurations), and radial and angular correlations in the 4f shell. Our calculated energy levels agree well with the experimental values. We have shown that low-lying highest-spin and second-highest spin states are reproduced very well, while for higher-lying states the accuracy of the calculation decreases. The procedure was verified by calculating optical emission spectra of NaYF₄:Eu,Tb; YAG:Eu,Tb; and Tb(acac)₃bpm (bpm is 2,2′-bipyridine, acac is acetylacetonate, and YAG is yttrium aluminium garnet). For these compounds ligand-field induced electric-dipole transition intensities were calculated.