  Density functional theory (DFT), including its extensions designed to treat
strongly correlated localized electron systems such as DFT+$U$ and
DFT+dynamical mean field theory, has proven exceedingly useful in studying the
magnetic properties of solids. However, materials with rare earths ($R$) have
remained a notable exception. The most vital rare-earth magnetic properties,
such as magnetocrystalline anisotropy (MA), have been notoriously elusive due
to the ubiquitous self-interaction error present in nearly all available DFT
flavors. In this work, we show explicitly how the orbital dependence of
self-interaction error may contradict Hund's rules and plague MA calculations,
and how analyzing DFT metastable states that respect Hund's rules can alleviate
the problem. We systematically investigate and discuss several
rare-earth-containing families, $R$Co$_5$, $R_2$Fe$_{14}$B, $R$Fe$_{12}$, and
$R$Mn$_6$Sn$_6$, to benchmark the MA calculations in DFT+$U$. For all compounds
we investigated, we found that our methodology reproduces the magnetic
easy-axis, easy-plane, and non-trivial easy-cone anisotropies in full agreement
with low-temperature experimental measurements. Besides the fully-numerical ab
initio approach, we further illustrate an efficient semi-analytical
perturbation method that treats the crystal field as a perturbation in the
limit of large spin-orbit coupling. This approach evaluates the rare-earth
anisotropy by assessing the dependence of crystal-field energy on
spin-quantization axis rotation using $4f$ crystal-field levels obtained from
non-spin-orbit calculations. Our analytical method provides a quantitative
microscopic understanding of the factors that control MA and can be used for
predicting new high-MA materials.
