Is there a complete description of the equations of motion of such a system? In particular, what would the Lagrangian and Hamiltonian for such a system be?
It depends on whether the particles are assumed to be points, or extended charged bodies.
If they are extended charged bodies, there is, as far as I know, no knowledge of a unique model of the system. The problem is that due to relativity the charged body cannot be idealized as a rigid body, but it is a system with infinity of degrees of freedom, like ball made of jelly. Mathematically complete description would require a model of motion and mutual internal forces between the charged parts of the particle. We do not have any convincing model of this. There are some published works that work with a more simplistic model, where the charged body is very regular ellipsoid that undergoes little or no deformation (Lorentz, Abraham, more recently Yaghijan and Medina are often cited) and are able to derive some conclusions about it, but all these calculations are of approximate character.
If the charged particles are points, point particles have only a handful of degrees of freedom and can be described by single position and velocity vectors. The situation is much simpler and this makes this kind of model much more attractive. There were papers by Fokker and Tetrode at the beginning of 20th century that show how a particular model of interacting particles, fully relativistic and in agreement with Maxwell's equations, may be formulated. Their formulation was focused on eliminating the EM fields from the description and used a variational principle to obtain the equations of motion of the particles directly from the action, without any middleman in form of EM field. However, this approach restricts the solutions to a highly special solutions of the Maxwell equations - so called half-retarded half-advanced fields.
More general formulation that does not require such restriction on the fields was first published, I think, by J. Frenkel in his paper
J. Frenkel, Zur Elektrodynamik punktfoermiger Elektronen, Zeits. f. Phys., 32, (1925), p. 518-534. http://dx.doi.org/10.1007/BF01331692
For a shorter, more easy-to-read account, see also
R. C. Stabler, A Possible Modification of Classical Electrodynamics, Physics Letters, 8, 3, (1964), p. 185-187. http://dx.doi.org/10.1016/S0031-9163(64)91989-4
It is true Frenkel too proposes half-retarded half-advanced solutions as particularly interesting since they allow for stable motion of hydrogen atom particles, but his formalism actually does not require them, it allows for any EM field that obeys Maxwell's equations.
The core idea is that particles act on other particles but never on themselves. The reason for this assumption for Frenkel was that self-action of a point on itself is contradictory and leads nowhere.
A particle acts on other particles via electromagnetic field of its own,so each field acquires an index that indicates which particle the field 'belongs to'. For example, particle $a$ generates electric field and its value at point $\mathbf r_b$ is $\mathbf E_a(\mathbf r_b)$. This is introduced so we can keep track of which field acts on which particle.
The fields obey the Maxwell equations with the owning particle as source:
$$
\nabla \cdot \mathbf E_a = \rho_a/\epsilon_0
$$
$$
\nabla \cdot \mathbf B_a = 0
$$
$$
\nabla \times \mathbf E_a = - \frac{\partial \mathbf B_a}{\partial t}
$$
$$
\nabla \times \mathbf B_a = \mu_0 \mathbf j_a + \mu_0\epsilon_0 \frac{\partial \mathbf E_a}{\partial t}
$$
Superposition of the elementary fields of all particles still obeys the Maxwell equations (thanks to their linearity), so this superposition is a good candidate for macroscopic total EM field.
The equation of motion of a charged particle $b$ is
$$
m_b \frac{d(\gamma_b \mathbf v_b)}{dt} = \sum_a' q_b \mathbf E_a(\mathbf r_b) + q_b\mathbf v_b \times \mathbf B_a(\mathbf r_b)
$$
(the prime near the sum sign means in case $a$ = $b$ the term is to be omitted)
This is a general formulation, fully relativistic and obeying both the Maxwell equations and the Lorentz force formula.
This direct formulation of equations of motion can be used to infer and check a variational Lagrangian formulation, where both field and particle variables are Lagrangian variables. The Lagrangian is
$$
L = \int d^3\mathbf x \mathcal{L}
$$
where
$$
\mathcal{L} = \sum_a\sum_b' -\frac{1}{4} F_a^{\mu\nu}F_{b,\mu\nu} + \sum_a\sum_b' j_a^\mu A_{b,\mu} - \sum_a m_a c^2\sqrt{1-v_a^2/c^2} \delta(\mathbf x - \mathbf r_a)
$$
Back to the case with two particles. One can make an approximation: the potentials are as if the particles were static or moving with speeds much lower than speed of light. These potentials can be inserted into the Lagrangian and then another Lagrangian, a function of particle positions and their derivatives, with no field variables, can be obtained.
In this way, the effect of fields is approximately expressed as a function of the particle variables. In the simplest case, this is the Coulombic term $\frac{q_aq_b}{4\pi\epsilon_0|\mathbf r_a - \mathbf r_b|}$. If the particle kinetic term of $L$ is linearized, one obtains non-relativistic Lagrangian of two particles interacting via static electric forces. For this approximate Lagrangian, one one can do the Legendre transformation and derive the common Hamiltonian function for the hydrogen atom. From the derivation it is clear that all this ignores magnetic interactions, retardation of interaction and is valid only for speeds much lower than the speed of light.
If better approximation is desired, one may insert potentials of particles moving rectilinearly with low speed and then the resulting Lagrangian function contains, in addition to the Coulombic term, a term that describes magnetic interaction. The term is called the Darwin interaction term and the whole Lagrangian the Darwin Lagrangian.
