  We apply the first-principles density functional theory + dynamical mean
field theory framework to evaluate the crystal field splitting on rare earth
sites in hard magnetic intermetallics. An atomic (Hubbard-I) approximation is
employed for local correlations on the rare earth 4$f$ shell and
self-consistency in the charge density is implemented. We reduce the density
functional theory self-interaction contribution to the crystal field splitting
by properly averaging the 4$f$ charge density before recalculating the
one-electron Kohn-Sham potential. Our approach is shown to reproduce the
experimental crystal field splitting in the prototypical rare earth hard magnet
SmCo$_5$. Applying it to $R$Fe$_{12}$ and $R$Fe$_{12}X$ hard magnets ($R=$Nd,
Sm and $X=$N, Li), we obtain in particular a large positive value of the
crystal field parameter $A_2^0\langle r^2\rangle$ in NdFe$_{12}$N resulting in
a strong out-of-plane anisotropy observed experimentally. The sign of
$A_2^0\langle r^2\rangle$ is predicted to be reversed by substituting N with
Li, leading to a strong out-of-plane anisotropy in SmFe$_{12}$Li. We discuss
the origin of this strong impact of N and Li interstitials on the crystal field
splitting on rare earth sites.
