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 noncollinear and high-rank intersite exchange, creating a complected configurational space with multiple minima. Although the $+\mathrm{U}$ correction to DFT allows, in principle, the modeling of such magnetic ground states, its results strongly depend on the initially symmetry breaking, constraining the nature of order parameter in the converged $\mathrm{DFT}+\mathrm{U}$ solution. As a rule, $\mathrm{DFT}+\mathrm{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 nonzero local contributions. Here we show how such phases can be efficiently captured using an educated constrained initialization of the on-site 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 $\mathrm{DFT}+\mathrm{U}$ methods. We assess the reliability of this procedure on the ferro-octupolar ground state recently predicted in ${\mathrm{Ba}}_{2}M{\mathrm{OsO}}_{6}\phantom{\rule{0.28em}{0ex}}(\mathrm{M}=\mathrm{Ca},\mathrm{Mg},\mathrm{Zn})$ [Phys. Rev. Lett. 127, 237201 (2021)].