Confusion between the Zeeman energy and electromagnetic energy The Zeeman term of the energy is something like
$$
\tag{1}
F_\text{Zeeman}=-\frac{1}{\mu_0}\int \vec H\cdot \vec M \,dV
$$
This term is in the Helmholtz free energy, I think; sadly I'm not sure. And the electromagnetic energy is something like
$$
\tag{2}
U=\frac12\int\left(\vec E\cdot \vec D + \vec H\cdot \vec B\right)\,dV
$$
This is deduced from the Lorentz force and Maxwell equations I think.
I'm confused because $\vec B=\mu_0 \vec H+\vec M$ and if we ignore the electric term we have
$$
U_\text{magn}=\frac12\int\vec H\cdot \left( \vec H+\tfrac{1}{\mu_0}\vec M\right)\,dV
=\frac12\int H^2\,dV +\frac{1}{2\mu_0}\int\vec H\cdot \vec M\,dV
$$
I'm confused because it seems to be almost $U_\text{magn}=-F_\text{Zeeman}$, so they basically have 2 different behavior. In $(1)$ we have a behavior where the magnetization tends to align with $\vec H$ when you minimize the energy, and in $(2)$ we have kind of the opposite. I'm sure there more than one thing where I'm wrong so any comments are welcome.
 A: First of all, your eq. (2) is not valid in general; it only applies for magnetic materials with linear constitutive relations ${\vec H}=\frac{1}{\mu}\vec{B}$ (or equivalently, $\vec{M}=\chi_{m}\vec{H}$).  Almost all everyday materials have linear electrostatic behavior, but this is not the case for the magnetic behavior, because of the existence of ferromagnets, for which $\vec{M}\neq0$ even when there is no external field.  However, that issue is separate from the qualitative issue that your question is mostly about.
The issue with the directions of the fields arises from the fact that, inside a magnetized material, the magnetization $\vec{M}$ and the fundamental magnetic field $\vec{B}$ point in opposite directions.  This is the opposite of what occurs electrically, but the behavior can already be seen in the difference between physical electric and magnetic dipoles.

While the fields of electric and magnetic dipoles have the same limiting forms far away, the near fields are totally different—and, in fact, point in opposite directions!  In the equatorial plane, the electric field of a physical dipole always points in the direction opposite the dipole moment $\vec{p}$.  However, the equatorial-plane magnetic field of a physical magnetic dipole (formed from a circulating current) points in the direction opposite the moment $\vec{m}$ outside the dipole by points parallel to $\vec{m}$ inside the dipole. (This is discussed more here:  Direction of electric dipole moment and magnetic dipole moment .)
For a magnetized material, built up of many tiny dipoles, this means that the field $\vec{B}$ points opposite the direction of magnetization $\vec{M}$ inside.

This means that the $-\vec{H}\cdot\vec{M}$ and $\vec{H}\cdot\vec{B}$ terms actually have the same signs inside a piece of (linearly) magnetized material.
A: Signs in formulae for the energy in the presence of a magnetic field, are often sources of headaches. Often the origin can be traced back to using the same name (energy) for quite different concepts and the same symbols for different fields.
Let me revisit your formulae to explain my previous statement.
Zeeman energy
For the moment, let's leave apart the free energy and let's write the term describing the interaction of a dipole moment $\vec M$ with an external magnetic field as
$$
\tag{1}
Z=\Delta E_\text{Zeeman}=-\frac{1}{\mu_0}\int \vec H_0\cdot \vec M \,dV=-\int \vec B_0\cdot \vec M \,dV.
$$
Notice, and this is an important point, that the fields $\vec B_0$ and $\vec H_0$ are the fields generated by the external sources in the absence of the dipole moment ($\vec H_0=\mu_0 \vec B_0 $). Moreover, in this formula, the $Z$ should be considered a function of $\vec B_0$. Therefore,
$$
\tag{2} {\mathrm d}Z=-\int {\mathrm d}\vec B_0\cdot \vec M \,dV.
$$
Equation $(1)$ can be obtained from equation $(2)$ by integrating over the external field $\vec B_0$ from $0$ to its final value.
Magnetic energy
From the Maxwell equations and the expression for the work in the presence of currents and a non-conservative electric field (Lorentz's force  does not enter in this story), it is easy to get an expression for the magnetic energy stored in the field configuration in the case of a linear relationship between the $\vec B$ and $\vec H$ fields:
$$
\tag{3}U_\text{magn}=\frac12\int\vec H\cdot \vec B\,dV.
$$
In the general case, it is better to use the more general expression (see Landau&Lifshitz Electrodynamics of Continuous Media)
$$
\tag{4} {\mathrm d} U_\text{magn}= \int\vec H\cdot {\mathrm d} \vec B\,dV.
$$
Again, Eqn. $(3)$ can be obtained from Eqn. $(4)$, in the case of a linear relationship between the two magnetic fields, by integrating over the field $\vec B$ from $0$ to its final value.
Notice that in formulae $(3)$ and $(4)$, the magnetic fields are not the fields in the absence of matter:
$$\tag{5}\vec B=\mu_0 (\vec H+\vec M)$$
Putting it all together
To compare Eqn.$(4)$ to Eqn.$(2)$, the first step is to re-express Eqn.$(4)$ in terms of the dipole moment $\vec M$ and the field $\vec B_0 $ (the external field in the absence of magnetic dipoles). This can be done (see Stratton's book) as follows:
$$
\begin{align}
\tag{6.1}\int\vec H\cdot {\mathrm d} \vec B\,dV&=\int\vec H_0\cdot {\mathrm d} \vec B_0\,dV+\\
\tag{6.2}&\int(\vec H_0\cdot {\mathrm d} \vec B- \vec B_0\cdot {\mathrm d} \vec H)\,dV+\\
\tag{6.3}&\int(\vec B_0\cdot {\mathrm d} \vec H_0- \vec H_0\cdot {\mathrm d} \vec B_0)\,dV+\\
\tag{6.4}&\int\vec B_0\cdot ({\mathrm d} \vec H-  {\mathrm d} \vec H_0)\,dV+\\
\tag{6.5}&\int(\vec H- \vec H_0)\cdot {\mathrm d} \vec B\,dV.
\end{align}
$$
The term $(6.3)$ vanishes because $\vec H_0=\mu_0 \vec B_0 $. Terms $(6.4)$ and $(6.5)$ are zero because integrals over the space of the scalar product of a solenoidal field (the $\vec B_0$ or the ${\mathrm d} \vec B$ fields) and irrotational fields ($\vec H_0$ and $\vec H$ have the same curl since originate from the same external sources). Finally, by using Eqn. $(5)$, the term $(6.2)$ can be recast as
$$
\tag{7}\int \vec B_0 \cdot {\mathrm d} \vec M\,dV.
$$
Eqn. $(7)$, provides an expression for the change of the energy-as a function (better as a functional) of $\vec M$. It also implies that
$$
\vec B_0({\bf r})=\frac{\delta U_\text{magn}}{\delta \vec M({\bf r})}.
$$
Therefore, we can go from a function of $\vec M$ to new energy depending on the external field $\vec B_0$ as independent variable, via a (functional) Legendre transform:
$$
\tilde U = U_\text{magn} - \int \vec B_0 \cdot \vec M
$$
so that the term corresponding to Eqn. $(7)$, when the external field is used as the independent variable becomes
$$
-\int {\mathrm d}\vec B_0\cdot \vec M \,dV,
$$
in agreement with Eqn. $(2)$.
To summarize, the two formulae are fully consistent. The reason for calling Eqn. $(2)$ a free energy is somewhat related to the fact that we have two different functions, related via a Legendre transform, to account for the energy variation of the system, according to which quantity we consider as the independent variable.
