Electromagnetic current for interaction with Dirac spinors The covariant form of the Dirac equation is given by
$$(i\gamma^{\mu}\partial_{\mu} - M) \Psi(x) = 0 $$
Einstein's summation is implied here, $x=(x^0,x^1,x^2,x^3)^T$.
I am simply looking for the Dirac theory equivalent to the classic electromagnetic interaction energy
$$ W = \int_V \rho \phi - \frac{1}{c} \mathbf{j\cdot A} \ dV $$
I have shown the Dirac equation's covariance and current conservation, but I can't seem to make the transition to the interaction energy.
Also, I have made myself familiar with the Klein-Gordon equation
$$(\square + M^2) \Psi(x) = 0$$
In that case, the interaction energy was "brought into" the equation by making the transition $ p^{\mu} \rightarrow p^{\mu} - \frac{e}{c} A^{\mu} $, where $p$ denotes the momentum operator, so the equation transforms to the form
$$(i\hbar\partial_{\mu} - \frac{e}{c}A_{\mu})(i\hbar\partial^{\mu} - \frac{e}{c}A^{\mu})\Psi(x) = M^2 \Psi(x)$$
Could it be that simple for the Dirac case, too?
 A: The equation of motion that you have is for a free Dirac field. That is, there is no photon (electromagnetic) field for the Dirac field to interact with. You would need to add an interaction term to your Lagrangian. 
Have you done interacting QFT? 
I believe adding a Maxwell Lagrangian 
$$ \mathcal{L}_M = −\frac{1}{4} \mathcal{F}_{\mu \nu} \mathcal{F}^{\mu \nu} $$
and promoting your partial derivative to a gauge covariant derivative 
$$ \partial_{\mu}  \rightarrow \mathcal{D}_{\mu} = \partial_{\mu} − i e A_{\mu}(x) $$
should give you what you want.
A: The probability current for the Dirac equation is $$j_p^\mu = \bar\psi \gamma^\mu \psi$$
and therefore the electromagnetic current is $$j^\mu =  -e\bar\psi\gamma^\mu\psi$$ since the charge of the Dirac field is $-e$. So to the Lagrangian we just add the usual term coupling current to vector potential $$\mathcal L = \bar\psi(i\gamma^\mu\partial_\mu - m)\psi - e\bar\psi\gamma^\mu\psi A_\mu.$$
Now the Euler-Lagrange equation for $\bar\psi$ can be read off as $$(i\gamma^\mu \partial_\mu - m -e \gamma^\mu A_\mu)\psi = 0.$$
(Different conventions will change some signs here but that's not really important.)
