Relativistic 2-body problem with instantaneous Colomb force Setup
Assume that we have two (isolated) massive and electrically charged particles, moving in the vacuum space, that interact through the instantaneous Coulomb force. I am an inertial observer and I want to build the equations of motion for each particle. Let ${{\mathbf{r}}_{i}}\left( t \right)$ be trajectory* of the particle with rest mass ${{m}_{i}}$ and electric charge ${{q}_{i}}$, and ${{\mathbf{F}}_{i}}\left( t \right)$ the force applied to that particle. Then, in my reference frame, the equations of motion for each particle would be:
$$\tfrac{d}{dt}{{\mathbf{p}}_{i}}\left( t \right)={{\mathbf{F}}_{i}}\left( t \right)$$
where ${{\mathbf{p}}_{i}}\left( t \right)\equiv {{m}_{i}}{{\gamma }_{i}}\left( t \right){{\mathbf{\dot{r}}}_{i}}\left( t \right)$, ${{\gamma }_{i}}\left( t \right)\equiv \gamma \left( {{v}_{i}}\left( t \right) \right)$ and ${{v}_{i}}\left( t \right)\equiv \left\| {{{\mathbf{\dot{r}}}}_{i}}\left( t \right) \right\|$.
${*}$ as it is observed in my reference frame, where $t$ is the reading of my clock
The question
Let, $\mathbf{r}\left( t \right)\equiv {{\mathbf{r}}_{1}}\left( t \right)-{{\mathbf{r}}_{2}}\left( t \right)$ and $K\equiv {{{q}_{1}}{{q}_{2}}}/{\left( 4\pi {{\varepsilon }_{0}} \right)}\;$. Which expression should I use for ${{\mathbf{F}}_{i}}\left( t \right)$? Possible candidates are:


*

*the Coulomb force between two particles that are ${r}\left( t \right)$ far apart from each other, i.e.:
$${{\mathbf{F}}_{1}}\left( t \right)=K\frac{1}{{{r}^{2}}\left( t \right)}\mathbf{\hat{r}}\left( t \right)=-{{\mathbf{F}}_{2}}\left( t \right)$$

*the Coulomb force that takes into account the phenomenon of length contraction, i.e.:
$${{\mathbf{F}}_{1}}\left( t \right)=K\frac{1}{{{\left( r\left( t \right){{\gamma }_{1}}\left( t \right) \right)}^{2}}}\mathbf{\hat{r}}\left( t \right)\quad ,\quad {{\mathbf{F}}_{2}}\left( t \right)=-K\frac{1}{{{\left( r\left( t \right){{\gamma }_{2}}\left( t \right) \right)}^{2}}}\mathbf{\hat{r}}\left( t \right)$$


If non of the above is true, which would be the correct expression? Let it be noted that we disregard any retardation effects.
My intuition
I would say with certainty that it is case 2., unless I had noticed the following:
$${{\mathbf{F}}_{1}}\left( t \right)+{{\mathbf{F}}_{2}}\left( t \right)=K\frac{1}{{{r}^{2}}\left( t \right)}\left( \frac{1}{\gamma _{1}^{2}\left( t \right)}-\frac{1}{\gamma _{2}^{2}\left( t \right)} \right)\mathbf{\hat{r}}\left( t \right)=K\frac{1}{{{c}^{2}}{{r}^{2}}\left( t \right)}\left( v_{2}^{2}\left( t \right)-v_{1}^{2}\left( t \right) \right)\mathbf{\hat{r}}\left( t \right)$$
which, in general, implies that $\tfrac{d}{dt}\left( {{\mathbf{p}}_{1}}\left( t \right)+{{\mathbf{p}}_{2}}\left( t \right) \right)\ne 0$, unless $v_{1}^{2}\left( t \right)=v_{2}^{2}\left( t \right)$ for any $t\in \mathbb{R}$. 
More questions
If my intuition is correct, then why, in general, is the total relativistic momentum not conserved? Is it because we neglected the retardation effects, and treated the problem semi-relativistically. Or is there a logical fallacy in the above considerations?
 A: Columb law (vetorial)
$$ \vec{E}(x,y,z) = \frac{e \vec{r}}{4 \pi \epsilon_0 r^3}, $$
Biot-Savart (vectorial + some calculus simplification from $ c = \frac{1}{\sqrt{\epsilon \mu}}$)
$$ \vec{B}(x,y,z) = \frac{\mu_0 e (\vec{v} \times \vec{r}) }{4 \pi r^3} =  
\frac{(\vec{v} \times \vec{r}) q}{c^2 (4 \pi \epsilon_0  r^3)} = \frac{\vec{v} \times \vec{E}}{c^2} $$
The correct way to treat electromagnetism interaction in relativity is to use
electromagnetic tensor $F^{\mu\nu}$
$$ F^\mu_\nu = 
\left( \begin{matrix}
0 & -E_x/c & -E_y/c & -E_z/c \\
-E_x/c & 0 & -B_z & B_y \\
-E_y/c & B_z & 0 & -B_x \\
-E_z/c & -B_y & B_x & 0 \\
\end{matrix} \right) ,
$$
Four-velocity vector
$$ \gamma = \frac{1}{\sqrt{1 - \frac{v^2}{c^2}}},  d\tau = \frac{dt}{\gamma}$$
$$ U^\mu = \left ( 
\begin{matrix} 
\gamma c \\ 
\gamma v_x \\ 
\gamma v_y \\ 
\gamma v_z \\
\end{matrix}\right ) , $$
Relativistic "force" law 
$$ \frac{d(m_{0}U^\mu)}{d\tau} = - e F^\mu_\nu U^\nu $$
Some notations
$$ \beta_x = \frac{v_x}{c} , \beta_x = \frac{v_y}{c} , \beta_z = \frac{v_z}{c} , \beta = \frac{v}{c}, v = \sqrt{v_x^2 + v_y^2 + v_z^2} $$
Lorentz Boost in arbitrary direction:
$$ \Lambda^\mu_\nu(v) = 
\left ( \begin{matrix}
\gamma    &  -\beta_x\gamma    &  -\beta_y\gamma   &   -\beta_z\gamma   \\
-\beta_x\gamma &   1 + (\gamma - 1)\beta_x^2/\beta^2  &   (\gamma - 1)\beta_x\beta_y/\beta^2   &      (\gamma - 1)\beta_x\beta_z/\beta^2 \\
-\beta_y\gamma &   (\gamma - 1)\beta_y\beta_z/\beta^2     & 1 + (\gamma - 1)\beta_y^2/\beta^2   &  (\gamma - 1)\beta_y\beta_z/\beta^2 \\
-\beta_z\gamma &    (\gamma - 1)\beta_z\beta_x/\beta^2    &   (\gamma - 1)\beta_z\beta_y/\beta^2    &   1 + (\gamma - 1)\beta_z^2/\beta^2 \\
\end{matrix} \right ) $$
Lowering + uppering indexes
$$ \Lambda_\mu^\nu(v) = \eta_{\mu\alpha}\eta^{\nu\beta} \Lambda^\alpha_\beta(v) $$
$$  \Lambda_\mu^\nu(v) = \left ( \begin{matrix}
\gamma    &  \beta_x\gamma    &  \beta_y\gamma   &   \beta_z\gamma   \\
\beta_x\gamma &   1 + (\gamma - 1)\beta_x^2/\beta^2  &   (\gamma - 1)\beta_x\beta_y/\beta^2   &      (\gamma - 1)\beta_x\beta_z/\beta^2 \\
\beta_y\gamma &   (\gamma - 1)\beta_y\beta_z/\beta^2     & 1 + (\gamma - 1)\beta_y^2/\beta^2   &  (\gamma - 1)\beta_y\beta_z/\beta^2 \\
\beta_z\gamma &    (\gamma - 1)\beta_z\beta_x/\beta^2    &   (\gamma - 1)\beta_z\beta_y/\beta^2    &   1 + (\gamma - 1)\beta_z^2/\beta^2 \\
\end{matrix} \right ) $$
Transformation of fields (we treat field as independent from source after generated - aka like a photon after being generated)
if $|\vec{v_{src}}| > 0$ then (src = source = "emitter", $v_{src}$ as seen by a rest observer)
$$ (F')^\mu_\nu = \Lambda^\mu_\alpha(v_{src}) \Lambda_\nu^\beta(v_{src}) F^\alpha_\beta $$ 
The full law:
$$ \frac{d(m_{0}U^\mu)}{d\tau} = - e (F')^\mu_\nu U^\nu $$

(In this image: the relativistic effect of precession - similar to the precession of the perihelion of Mercury [as opposed to classical prediction which would've been a simple unique orbit])
PS: If you need to create a simulation of some sort I recommend you to use matrix(vectorial) calculus and not try to decompose(break/dismantle) tensors into simpler forms because it's for sure you'll end up making a sign/symbol mistake. 
References (some):


*

*https://en.wikipedia.org/wiki/Covariant_formulation_of_classical_electromagnetism 

*A Boost in an Arbitrary Direction

*The Feynman Lectures on Physics, Volume II - Lorentz Transformations of the Fields

*Nicolae Barbulescu - Bazele Fizice ale Relativitatii Einsteiniene, 1975

*A.Einstein - The Meaning of Relativity, Sixth Edition

