  Using the multipolar expansion of electrostatic and magnetostatic potential
energies, we characterize the long-range interactions between two weakly-bound
diatomic molecules, taking as an example the paramagnetic Er$_2$ Feshbach
molecules which were produced recently. Since inside each molecule, individual
atoms conserve their identity, the intermolecular potential energy can be
expanded as the sum of pairwise atomic potential energies. In the case of
Er$_2$ Feshbach molecules, we show that the interaction between atomic magnetic
dipoles gives rise to the usual $R^{-3}$ term of the multipolar expansion, with
$R$ the intermolecular distance, but also to additional terms scaling as
$R^{-5}$, $R^{-7}$, and so on. Those terms are due to the interaction between
effective molecular multipole moments, and are strongly anisotropic with
respect to the orientation of the molecules. Similarly the atomic pairwise van
der Waals interaction results in $R^{-6}$, $R^{-8}$, ... terms in the
intermolecular potential energy. By calculating the reduced electric-quadrupole
moment of erbium ground level $\langle J=6||\hat{Q}_2||J=6\rangle = -1.305$
a.u., we also demonstrate that the electric-quadrupole interaction energy is
negligible with respect to the magnetic-dipole and van der Waals interaction
energies. The general formalism presented in this article can be applied to
calculate the long-range potential energy between arbitrary charge
distributions composed of almost free subsystems.
