  Multipolar magnetic phases in correlated insulators represent a great
challenge for Density Functional Theory (DFT) due to the coexistence of
intermingled interactions, typically spin-orbit coupling, crystal field and
complex non-collinear and high-rank inter-site exchange, creating a complected
configurational space with multiple minima. Though the +U correction to DFT
allows, in principle, the modelling of such magnetic ground states, its results
strongly depend on the initially symmetry breaking, constraining the nature of
order parameter in the converged DFT+U solution. As a rule, DFT+U calculations
starting from a set of initial on-site magnetic moments result in a
conventional dipolar order. A more sophisticated approach is clearly needed in
the case of magnetic multipolar ordering, which is revealed by a null integral
of the magnetization density over spheres centered on magnetic atoms, but with
non-zero local contributions. Here we show how such phases can be efficiently
captured using an educated constrained initialisation of the onsite density
matrix, which is derived from the multipolar-ordered ground state of an ab
initio effective Hamiltonian. Various properties of such exotic ground states,
like their one-electron spectra, become therefore accessible by all-electron
DFT+U methods. We assess the reliability of this procedure on the
Ferro-Octupolar ground state recently predicted in Ba$_2$MOsO$_6$ (M = Ca, Mg,
Zn) [Phys. Rev. Lett. 127, 237201 (2021)]
