Is it possible to propagate a relativistic system of particles in time using Verlet? The Verlet algorithm and its derivations are very popular methods to integrate Newton's equations of motion in time and obtain a trajectory for a system with $N$ particles. 
I work with classical molecular mechanics in several of my research projects, but they all apply force fields as the potential that govern the motion of the $N$ bodies. Recently, however, I was thinking if the Verlet Algorithm (with which I have some familiarity) could be used to integrate some other equations of motion and to obtain trajectory which take into account effects of special relativity such as time dilation and length contraction.
I have a good grasp of overall relativity but I am a theoretical chemist, not a physicist and my usual realm of practise is quantum chemistry, so it turns out I don't have enough mathematical expertise and/or practise with relativity numerical idiossyncrasies to figure that out on my own, so I came to seek wisdom. 
So my question regarding all that is threefold:
$(1)$ What equations of motion should be integrated to yield the trajectories?
$(2)$ Can Verlet algorithm be used to integrate them?
$(3)$ If not, then which method should I use?
Given that the $N$ particles are subject to pair-wise inverse-square potentials.
 A: Verlet implicitly assumes forces are a function of position only. This means verlet is not best suited to forces that are a function of velocity as well as position. For example, drag. Verlet is not the best of choices in cases where drag is significant. The second order Newmark-beta method works better than verlet in cases where drag is significant.
The same applies to relativistic condition because relativistic force involves velocity. The problem with Newmark-beta is that it is an implicit method. You can't magically determine what the forces are at the end of the time interval. You can only guess, and then iterate until things converge. That's what the Newmark-beta method do. One iteration might be enough, but even then there's an extra step over verlet.
Side comment: Verlet and Newmark-beta are not the be-all and end-all of integration techniques. You use verlet because you are a chemist. Your N-body problems inherently involve a rather large N. You have but no choice to use a low order technique. Anything else would kill you computationally. Low order techniques such as verlet (and Newmark-beta) are rather suboptimal when the number of bodies is small.
Update
Just to be clear, I'm addressing this question from the perspective of various astrophysics / astronomical / aerospace engineering integrators for systems with a small number (three to hundreds) of bodies. Which technique one uses depends on a number of factors. Just a few are:


*

*The time span of the integration.
The conservative nature of symplectic integrators is important if the time span is very long (millions of years). People are continuously improving exiting techniques and developing new techniques. Here's a link to a scholar.google.com search on symplectic partitioned Runge-Kutta (SPRK) integrators published just during this year (2014). Symplecticity is not important at all if the time span is very short, where accuracy takes precedence.

*The desired accuracy.
Up to a point, higher order techniques tend to offer much greater accuracy. Most techniques have a problem as order increases. Increasing the order results in ever more and ever larger terms with oscillating signs. For example, using 16th order Gauss-Jackson with double precision arithmetic makes no sense whatsoever. It does make sense if one uses quad precision, but now there's a huge performance hit because most floating point processors don't implement true quad precision. (Gauss-Jackson is essentially an Adam-style technique specialized for second order ODEs.) JPL uses some other Adams-style technique to develop it's Development Ephemeris.

*Whether collisions occur.
Collisions (or rather, close interactions) bedevil high order techniques. The purpose of using high order techniques is that this enables huge time steps. In fact, you want to take huge time steps with a high order technique. Up to a point, bigger steps means more accuracy. Ideally you want the integration to live just at the cusp where errors due to use of finite precision arithmetic and errors inherent to the integration technique (truncation errors) are in balance.

*Whether there's a large central body.
There are a large number of specialized techniques that apply when there's a single solar mass gorilla in the room. Some of these are very old and predate computers by hundreds of years. Lagrange's planetary equations (and many variants, e.g., Gauss' planetary equations) are obviously old, but still useful. They yield time derivatives of Keplerian elements given small perturbations on top of the gravitation by that solar mass gorilla. People still use these techniques. Another oldie-but-goodie is Encke's method. That solar mass gorilla makes objects follow a perturbed Keplerian orbit. Encke-type methods follow a reference Keplerian segment and compute perturbations on it. The reference trajectory is regularly readjusted when the perturbed trajectory deviates to far from the reference trajectory.

*How many computers one has on hand.
There are a number of ways to partition the problem across multiple computers. Exchange states, compute accelerations, exchange accelerations, take a step. With a small N, this can be done with coarse grain parallelism. With a moderately large N, this can be done using a computer's graphics card. Much better to put that graphic card to work doing useful computations than running a computer game!

*Constraints on step size.
Spacecraft turn their thrusters on and off with alarming frequency. Many flight software systems have a multi-hertz cycle rate. A spacecraft simulation with flight software running in the loop means the step size is dictated by that cycle rate. There's no point in using a high order technique that likes to take minute-long or longer steps if the flight software turns jets on and off at 100 hertz.
A: (This answer really only gives the relevant ODEs).
To answer 1., I'll suppose you're given a potential. This link derives the relativistic Lagrangian $L=-\sqrt{1-\beta^2}-U$ giving $\frac{d}{dt}\left( \gamma \vec{\beta}\right)=-\nabla U$, where I've set $c=1$, $m=1$, and $\gamma=\frac{1}{\sqrt{1-\beta^2}}$. If $U$ reacts in some way to the motion of the particle (which it must in reality, since particles can't interact over a distance) I believe that you'd have to use something like retarded Lienerd-Wiechert potentials.
I don't have an answer to 2. In one dimension you can do this:


*

*Write $\beta=\tanh(k)$

*Then, using the identity $1-\tanh^2(k)=\mbox{sech} ^2(k)$, $\gamma \beta=\sinh(k)$

*The formula is them $\dot{k}=-\frac{\nabla U}{\cosh(k)}$


You might be able to do something similar in 3D, but in the end we're after positions and it may or may not be the case that switching between rapidity and velocity is worth it.
