The equation of motion is given by (in Heaviside-Lorentz units)
[tex]m \frac{\mathrm{d}}{\mathrm{d} t} \left (\frac{\dot{\vec{x}}}{\sqrt{1-\dot{\vec{x}}/c^2}} \right)=q \left (\vec{E}+ \frac{\vec{v}}{c} \times \vec{B} \right )=q \left [-\vec{\nabla} \phi + \frac{\vec{v}}{c} \times (\vec{\nabla} \times \vec{A}) \right ].[/tex]
To derive it from the principle of least action, I prefer the Lagrangian form, because there the action is manifestly covariant:
[tex]L=-mc^2 \sqrt{1-\vec{v}^2/c^2}-q \left (\phi-\frac{\dot{\vec{x}}}{c} \cdot \vec{A} \right ).[/tex]
The action is manifestly covariant, because
[tex]\sqrt{1-v^2/c^2} \mathrm{d} t = \mathrm{d} \tau = \sqrt{\mathrm{d} x_{\mu} \mathrm{d} x^{\mu}}[/tex]
and
[tex](\phi-\frac{\dot{\vec{x}}}{c} \cdot \vec{A}) \mathrm{d} t = \mathrm{d} t \frac{1}{c} \frac{\mathrm{d} x^{\mu}}{\mathrm{d} t} A_{\mu}=\frac{1}{c} \mathrm{d} x^{\mu} A_{\mu}.[/tex]
The four potential [itex](A^{\mu})=(\phi,\vec{A})[/itex] is a Mikowski four-vector.
To get the equations of motion, you just write down the Euler-Lagrange equations
[tex]\frac{\mathrm{d}}{\mathrm{d} t} \frac{\partial L}{\partial \dot{\vec{x}}}=\frac{\partial L}{\partial \vec{x}}.[/tex]
Just taking the derivatives, gives
[tex]\vec{P}=\frac{\partial L}{\partial \dot{\vec{x}}}=\frac{m \dot{\vec{x}}}{\sqrt{1-\dot{\vec{x}}^2/c^2}}+\frac{q}{c} \vec{A}.[/tex]
Here it is important to note that [itex]\vec{P}[/itex] is the canonical momentum, not the mechanical momentum, which is given by the first term only.
Further we have
[tex]\frac{\partial L}{\partial \vec{x}}=-q \vec{\nabla} \phi + \frac{q}{c} \vec{\nabla} (\dot{\vec{x}} \cdot \vec{A}).[/tex]
Written in the three-dimensional Ricci calculus (with all indices as subscripts) the equation of motion thus reads
[tex]\dot{P}_j=\frac{\mathrm{d}}{\mathrm{d} t} \left (\frac{m \dot{x}_j}{\sqrt{1-\dot{\vec{x}}^2/c^2}} + q A_j \right ) = \frac{\mathrm{d}}{\mathrm{d} t} \left (\frac{m \dot{x}_j}{\sqrt{1-\dot{\vec{x}}^2/c^2}} \right ) + \frac{q}{c}\dot{x}_k \partial_k A_j = -q \partial_j \phi + \frac{q}{c} \dot{x}_k \partial_j A_k.[/tex]
Combining the two terms with the vector potential on the right-hand side gives
[tex]\frac{\mathrm{d}}{\mathrm{d} t} \left (\frac{m \dot{x}_j}{\sqrt{1-\dot{\vec{x}}^2/c^2}} \right )=-q \partial_j \phi + \frac{q}{c} \dot{x}_k (\partial_j A_k-\partial_k A_j).[/tex]
The last term can be rewritten as
[tex]\dot{x}_k (\partial_j A_k-\partial_k A_j) =\dot{x}_k \epsilon_{jkl} (\vec{\nabla} \times \vec{A})_l=[\dot{\vec{x}} \times (\vec{\nabla} \times \vec{A})]_j.[/tex]
Rewriting everything in vector notation yields the correct equation of motion.
Note that the notation in Wikipedia, from where you copied your equation, is not very clear. That's why I used the Ricci calculus in the intermediate steps, where it is clear what are free indices and which are repeated indices to be summed.