  In insulating lanthanides, unquenched orbital momentum and weak crystal-field
(CF) splitting of the atomic $J$ multiplet at lanthanide ions result in a
highly ranked (multipolar) exchange interaction between them and a complex
low-temperature magnetic order not fully uncovered by experiment. Explicitly
correlated {\it ab initio } methods proved to be highly efficient for an
accurate description of CF multiplets and magnetism of individual lanthanide
ions in such materials. Here we extend this {\it ab initio } methodology and
develop a first-principles microscopic theory of multipolar exchange
interaction between $J$-multiplets in $f$ metal compounds. The key point of the
approach is a complete account of Goodenough's exchange mechanism along with
traditional Anderson's superexchange and other contributions, the former being
dominant in many lanthanide materials. Application of this methodology to the
description of the ground-state order in the neodymium nitride with rocksalt
structure reveals the multipolar nature of its ferromagnetic order. We found
that the primary and secondary order parameters (of $T_{1u}$ and $E_g$
symmetry, respectively) contain non-negligible $J$-tensorial contributions up
to the ninth order. The calculated spin-wave dispersion and magnetic and
thermodynamic properties show that they cannot be simulated quantitatively by
confining to the ground CF multiplet on the Nd sites. Our results demonstrate
that the {\it ab initio } approach to the low-energy Hamiltonian represents a
powerful tool for the study of materials with complex magnetic order.
