I thought about your problem for a bit andd have arrived at some calculation that I would like to show you.
Note that we are using natural units where $c=1$.
Steps
- We will first determine the exact representation of $T^{EM}_{\mu\nu}$.
- We will then represent the formulation of the field equations of both the metric field and the elctromagnetic field.
- At last, we will derive the equation(s) of motion for a particle in these fields.
Finding the value of $T^{EM}_{\mu\nu}$
$T_{\mu\nu}$ can be represented in the following way:
$$T_{\mu\nu} = \frac {-2}{\sqrt{-g}} \frac {\delta\sqrt{-g}\mathscr{L}_{matter}}{\delta g^{\mu\nu}}$$
where $\mathscr{L}_{matter}$ is the Lagrangian density of the "matter field". This matter field is actually any kind of Lagrangian of any kind of field. Because we are working with EM fields, we will have to use:
$$\mathscr{L}_{matter} = -\frac{1}{4}F^{\mu\nu}F_{\mu\nu} + J^{\mu}A_{\nu}.$$
Note the use of Einstien's summation convention in this and the following equations.
The calculations are left to the reader, but after we have done all the awful calculations, we arrive at
$$T_{\mu\nu}^{EM} = -g^{\gamma\alpha}F_{\gamma\alpha}F_{\mu\nu} - \frac{1}{4}g_{\mu\nu} + g_{\mu\nu}J^{\alpha}A_{\alpha}$$
which is quite a clean expression.
Field equations
We have written down the exact formulation of $T_{\mu\nu}^{EM}$. The following steps are very easy.
We just write $T_{\mu\nu} = T_{\mu\nu}^{other} + T_{\mu\nu}^{EM}$ and we are good to go. Just apply this formulation of $T_{\mu\nu}$ to the Einstien's field equations and we could derive some solutions.
For the elctromagnetic field, we know the field equations look like this:
$$\partial_{\mu}F^{\mu\nu} = J^{\mu}.
$$
Surprisingly, to include the metric tensor, we just have to slightly tweak this equation(s) and use the covariant derivative instead of the normal derivative and we have
$$\nabla_{\mu}F^{\mu\nu} = J^{\mu}.$$
Hence, we have determined all the field equations.
The equation(s) of motion
The Lagrangian of a particle in a metric field and electromagnetic field are as follows:
$$L_{metric} = -mg_{\mu\nu}\frac{\mathrm{d}X^{\mu}}{\mathrm{d}t}\frac{\mathrm{d}X^{\nu}}{\mathrm{d}t} \\ L_{EM} = -m\sqrt{1 - \dot{X}^{m}\dot{X}_{m}} + \frac{\mathrm{d}X^{\mu}}{\mathrm{d}t}A_{\mu}.$$
Note that the greek indices on $X$ means that we are summing over all the dimensions including the time and the latin index means that we are summing only over the spatial dimensions.
The reader has to apply the Euler-Lagrange equation to $L_{metric} + L_{EM}$ and he will get the equation(s) of motion as described in the question:
$$\ddot{X}^{\rho} + \Gamma^{\rho}_{\mu\nu}\dot{X}^{\mu}\dot{X}^{\nu} = \frac{q}{m}g^{\rho\mu}F_{\mu\nu}\dot{X}^{\nu}.$$
To conclude, I am willing to listen to the clarifications of the readers and hope this helps.
Edit
Thanks for the clarification of the expert reviewers in the comments, I would like to mention some of the problems in my answer.
The $T^{EM}_{\mu\nu}$ is gauge variant, which means it can violate some gauge symmetries of the universe.
The last two lagrangians are for a point particle, no assumption should be made that it is true for any kind of matter distribution. The reader can change the lagrangian based on the problem they are solving.