A class of explicit numerical schemes is developed to solve for the relativistic dynamics and spin of particles in electromagnetic fields, using the Lorentz-Bargmann-Michel-Telegdi equation formulated in the Clifford algebra representation of Baylis. It is demonstrated that these numerical methods, reminiscent of the leapfrog and Verlet methods, share a number of important properties: they are energy conserving, volume conserving, and second-order convergent. These properties are analyzed empirically by benchmarking against known analytical solutions in constant uniform electrodynamic fields. It is demonstrated that the numerical error in a constant magnetic field remains bounded for long-time simulations in contrast to the Boris pusher, whose angular error increases linearly with time. Finally, the intricate spin dynamics of a particle is investigated in a plane-wave field configuration.