In Dirac's Lectures on Quantum Mechanics he writes down the Hamiltonian for the free electromagnetic field:
\begin{equation} H = \int \frac{1}{4} F^{rs} F_{rs} + \frac{1}{2} \pi^r \pi^r - A_0 \pi^r{}_{,r} \, \mathrm{d}^3 x \end{equation}
(Note that Dirac uses $B$ to denote the conjugate momentum, so I have replaced it with $\pi$, which is less confusing.)
If one evaluates the first two terms, one recovers the familiar expression for the energy density of the field, namely $\frac{1}{2}(E^2 + B^2)$. The third term will eventually be seen to vanish once Gauss's law is recovered.
My question is: if we started out with the full Lagrangian which includes the interaction term $-\frac{1}{c} J^\mu A_\mu$, then what happens to this term when the Legendre transformation is performed? As far as I can tell, this term simply ends up in the result with changed sign, i.e.,
\begin{equation} H = H_m + \int \frac{1}{4} F^{rs} F_{rs} + \frac{1}{2} \pi^r \pi^r + \frac{1}{c} J^\mu A_\mu - A_0 \pi^r{}_{,r} \, \mathrm{d}^3 x \end{equation}
where $H_m$ denotes the Hamiltonian for matter particles alone. This now seems to give the wrong value of the energy density; the last term will cancel the $\rho \varphi$ term in $\frac{1}{c} J^\mu A_\mu$, but the $-\frac{1}{c} \mathbf{A} \cdot\mathbf{J}$ part will be left uncancelled. That means the Hamiltonian density will differ by $-\frac{1}{c} \mathbf{A} \cdot \mathbf{J}$ from the energy density
\begin{equation} u = u_m + \frac{1}{2}(E^2 + B^2) \end{equation}
So what part of this analysis is incorrect?
(Apologies if I got the units wrong; I am not proficient in cgs units.)