We know that the EM-Lagrangian \begin{align} -\frac{1}{16 \pi} F_{\mu \nu} F^{\mu \nu} - \frac{1}{c} j_{\mu} A^{\mu} \end{align} correctly describes the evolution of the electromagnetic field in the presence of a given $j(x) = (c \rho(\vec{x}, t), \vec{j}(\vec{x}, t))$, which is externally specified, and not governed by the Lorentz-force (or rather its equivalent action on charge densities). Treating $j^{\mu}$ as an independent field won't work, because then the e.l. equations would yield $A^{\mu} = 0$.
Is there any classical (classical as in: not quantized, no operators) lagrangian that yields the said equations, not only for the E.M. field, but also for the "sources"?
The next best thing I can think of is the EM-field coupled to a scalar field, with mass- and kinetic term for both and $j^{\mu} = \phi^* \partial_{\mu} \phi + \phi \partial^{\mu} \phi^*$ in the interaction term: \begin{align} -\frac{1}{16 \pi} F_{\mu \nu} F^{\mu \nu} - \frac{1}{c}( \phi^* \partial_{\mu} \phi + \phi \partial^{\mu} \phi^*) A_{\mu} + (\partial_{\mu} - A_{\mu})\phi^*(\partial^{\mu} - A^{\mu})\phi + (m^2)\phi^* \phi \end{align}
But is there a lagrangian that is even closer? If not, why isn't this classical lagrangian used to do classical calculations? Doesn't it reproduce the dynamics that would follow from the Lorentz force?
To make more explicit what the lagrangian should yield: It should yield Maxwell's equations, coupled to sources, and additionally the behaviour consistent with a charge density that describes point particles with a trajectory $r_i$ by $j_i^{\mu}(x) = (c, \vec{v}_i(t))^T \delta(\vec{x} - \vec{x}_i(t))$.