The covariant formalism is a bit tricky. The problem is that you cannot use the proper time directly in the Lagrangian, because it implicitly contains the constraint that \mathrm{d} \tau=\mathrm{d} t \sqrt{1-v^2} (I'm setting the speed of light to 1 here).
You can, however use an arbitrary "world parameter" and write down a parameter-independent action functional. For the free particle it's
S_0[x]=-m \int \mathrm{d} \lambda \sqrt{\dot{x}_{\mu} \dot{x}^{\mu}}.
I'm using the "west-coast convention" for the pseudometric, i.e., \eta_{\mu \nu}=\text{diag}(1,-1,-1,-1).
The interaction with an external electromagnetic field is determined by
S_i[x]=-q \int \mathrm{d} \lambda \dot{x}^{\mu} A_{\mu}(x).
Note that the total Lagrangian,
L=-m \sqrt{\dot{x}_{\mu} \dot{x}^{\mu}} - q \dot{x}^{\mu} A_{\mu}(x)
is a homogeneous function of degree 1 wrt. \dot{x}. This implies the parameter independence. Now you can derive the equations of motion from the Hamilton principle of least action as usual, using the Euler-Lagrange equations. The canonical momenta are given by
p_{\mu}=\frac{\partial L}{\partial \dot{x}^{\mu}}=-m \frac{\dot{x}_{\mu}}{\sqrt{\dot{x}_{\nu} \dot{x}^{\nu}}}-q A_{\mu}
and thus the equation of motion
\dot{p}_{\mu}=\frac{\partial L}{\partial x^{\mu}} \; \Rightarrow \; -m \frac{\mathrm{d}}{\mathrm{d} \lambda} \left (\frac{\dot{x}_{\mu}}{\sqrt{\dot{x}_{\nu} \dot{x}^{\nu}}} \right ) - q \dot{x}^{\nu} \partial_{\nu} A_{\mu}(x)=-q \dot{x}^{\nu} \partial_{\mu} A_{\nu}.
Now you can bring this into a more familiar form by rearranging the terms to
m \frac{\mathrm{d}}{\mathrm{d} \lambda} \left (\frac{\dot{x}_{\mu}}{\sqrt{\dot{x}_{\nu} \dot{x}^{\nu}}} \right ) = q \dot{x}^{\nu} (\partial_{\mu} A_{\nu}-\partial_{\nu} A_{\mu})=q F_{\mu \nu} \dot{x}^{\nu}.
Now you can choose for \lambda whatever parameter you like. It is crucial to note that the equation of motion, if derived from a Lagrangian that is homogeneous in \dot{x}^{\mu} of degree one, automatically fulfills the constraint equation
\dot{x}^{\mu} \frac{\mathrm{d}}{\mathrm{d} \lambda} \left (\frac{\dot{x}_{\mu}}{\sqrt{\dot{x}_{\nu} \dot{x}^{\nu}}} \right )=0,
which implies that the 4 equations of motion are not independent from each other and thus you can choose the parameter \lambda as you like after the variation is done, i.e., on the level of the equations of motion.
If you choose the proper time, \lambda=\tau, then you get \dot{x}_{\mu} \dot{x}^{\mu}=1 and thus
m \frac{\mathrm{d} u_{\mu}}{\mathrm{d} \tau}=q F_{\mu \nu} u^{\nu} \quad \text{with} \quad u^{\mu}=m \frac{\mathrm{d} x^{\mu}}{\mathrm{d} \tau}.
You can as well chose a coordinate time, t of an inertial frame. Then the equation of motion appears in a form that is not longer manifestly covariant, but you get, because of \dot{x}^{\mu} \dot{x}_{\nu}=1-\vec{v}^2, \quad \vec{v}=\frac{\mathrm{d} \vec{x}}{\mathrm{d} t},
for the spatial part of the equations of motion (writing the Faraday tensor in terms of the three-dimensional notation with \vec{E} and \vec{B})
m \frac{\mathrm{d}}{\mathrm{d} t} \left (\frac{\vec{v}}{\sqrt{1-\vec{v}^2}} \right ) = q (\vec{E}+\vec{v} \times \vec{B}),
and the time component is just the energy-work relation, following from this:
m \frac{\mathrm{d}}{\mathrm{d} t} \left (\frac{1}{\sqrt{1-\vec{v}^2}} \right ) = \vec{v} \cdot \vec{E}.
