What follows is what you could describe as naiively applying the Legendre transform for constructing the Hamiltonian from the Lagrangian. Weinberg's "Quantum Theory of Fields" Vol. I chapter 8 goes over the canonical quantization of the electromagnetic field that correctly handles constraint equations.
Rather than finding the $T^{00}$ component of the stress-energy tensor, which often needs to be modified using the equations of motion to bring them into a recognizable form, you can construct the Hamiltonian using the canonical formalism. First, take the Lagrangian (using that instead of the density so that I can be free to move back and forth from mode space - also using Heaviside-Lorentz units with $\hbar=c=1$, and particle physicists' metric signature $(+,-,-,-)$),
\begin{align}
L & = \int \operatorname{d}^3x \left(\bar{\psi}\left[i\gamma^\mu D_\mu - m\right]\psi - \frac{1}{4}F_{\mu\nu}F^{\mu\nu} \right) \\
& = \int \operatorname{d}^3x \left( \bar{\psi}\left[i\gamma^\mu D_\mu - m\right]\psi + \frac{1}{2} E_iE_i - B_iB_i \right) \\
& = \int \operatorname{d}^3x \left( \vphantom{\frac{1}{2}} \bar{\psi}i\gamma^0\partial_0\psi + \bar\psi i \gamma^i D_i \psi - m\bar{\psi}\psi - e \phi \bar{\psi}\gamma^0\psi \right. \\
& \hphantom{=\int \operatorname{d}^3x}\ \ \left. + \frac{1}{2}\left[-\partial_i\phi-\partial_0A_i\right] \left[-\partial_i\phi-\partial_0A_i\right] - \frac{\epsilon_{ijk}\epsilon_{inm}}{2} \partial_jA_k \partial_nA_m \right).
\end{align}
The reason for expanding the the covariant inner products is to make explicit the next step: defining the canonically conjugate momenta.
\begin{align}
\pi_\psi &\equiv \frac{\delta L}{\delta \partial_0 \psi} = \bar{\psi}i\gamma^0 = i\psi^\dagger \\
\Pi_i &\equiv \frac{\delta L}{\delta \partial_0 A_i} = \partial_i\phi + \partial_0 A_i \\
\Pi_\phi &\equiv \frac{\delta L}{\delta \partial_0 A_0} = 0
\end{align}
Notice how the momentum canonically conjugate to $\phi=A_0$ doesn't appear in the Lagrangian. This means that, without gauge fixing, $\phi$ is a non-dynamical field that doesn't have any momentum. In classical physics, it plays the roll of a Lagrange multiplier that enforces $\rho - \nabla\cdot\mathbf{E}= e\bar{\psi}\gamma^0\psi + \partial_i \Pi_i =0$, the first of Maxwell's equation.
Now, $H$ is defined as
\begin{align}
H &= \int \operatorname{d}^3x \left[\pi_\psi \partial_0\psi + \Pi_i \partial_0 A_i\right] - L \\
&= \int \operatorname{d}^3x \left[ \vphantom{\frac{\epsilon_i}{2}} \Pi_i (\Pi_i - \partial_i\phi) - \pi_\psi \gamma^0 \gamma^i D_i\psi -i m\pi_\psi \gamma^0\psi - ie\phi\pi_\psi \psi \right. \\
&\hphantom{= \int \operatorname{d}^3x}\ \ \left.- \frac{1}{2} \Pi_i \Pi_i + \frac{\epsilon_{ijk}\epsilon_{inm}}{2}\partial_jA_k \partial_nA_m \right].
\end{align}
Notice that the term containing the time derivative of $\psi$ vanishes entirely because the Dirac equation is first order in time, so all of the needed time derivatives are supplied by the canonical equations of motion.
Now we move the gauge invariant parts of the electromagnetic field into mode space (Fourier transform over space, but not time) to highlight some interesting structure in the Hamiltonian.
\begin{align}
H & = \int \operatorname{d}^3k \left[ \frac{1}{2}\Pi_i\Pi_i + \frac{k^2}{2} \left(\delta_{ij} - \frac{k_i k_j}{k^2}\right)A_iA_j \right] \\
&\hphantom{=}+\int \operatorname{d}^3x \left[- \pi_\psi \gamma^0 \gamma^i D_i\psi -i m\pi_\psi \gamma^0\psi - ie\phi\pi_\psi \psi - \Pi_i\partial_i\phi \right] \\
&= \int \operatorname{d}^3k \left[ \frac{k_ik_j}{2 k^2}\Pi_i\Pi_j + \frac{1}{2} \left(\delta_{ij} - \frac{k_i k_j}{k^2}\right)\left(\Pi_i\Pi_j + k^2 A_iA_j\right) \right] \\
&\hphantom{=}+\int \operatorname{d}^3x \left[- \pi_\psi \gamma^0 \gamma^i D_i\psi -i m\pi_\psi \gamma^0\psi - ie\phi\pi_\psi \psi + \phi \partial_i \Pi_i\right]
\end{align}
Note that the electromagnetic field momentum, $\Pi_i$, is gauge invariant (it's just $-E_i$). The quantity $\left(\delta_{ij} - \frac{k_i k_j}{k^2}\right)A_j$ is also gauge invariant since it is the Fourier transform of the solenoidal part of $A_i$; in other words, the purpose of that last line was to collect the solenoidal components of the electromagnetic field (the photons) and move the derivative from $\phi$ to $\Pi_i$ in that term. The longitudinal kinetic term, $\frac{k_i k_j}{2k^2}\Pi_i\Pi_j$, is likewise nicely isolated.
It is also interesting to note that the only terms that aren't manifestly gauge invariant are
$$-ie\phi\pi_\psi \psi + \phi \partial_i \Pi_i = e\phi\bar{\psi}\gamma^0\psi - \phi\partial_i E_i,$$
which vanishes by the equations of motion. It strikes me as likely that it is dealing with these terms in a canonical fashion that leads to the requirement for gauge fixing and constraints.
