Direct observations of compact objects, in the form of radiation spectra, gravitational waves from VIRGO/LIGO, and forthcoming direct imaging, are currently one the primary source of information on the physics of plasmas in extreme astrophysical environments. The modeling of such physical phenomena requires numerical methods that allow for the simulation of microscopic plasma dynamics in presence of both strong gravity and electromagnetic fields. In \cite{bacchini2018a} we presented a detailed study on numerical techniques for the integration of free geodesic motion. Here we extend the study by introducing electromagnetic forces in the simulation of charged particles in curved spacetimes. We extend the Hamiltonian energy-conserving method presented in \cite{bacchini2018a} to include the Lorentz force and we test its performance compared to that of standard explicit Runge-Kutta and implicit midpoint rule schemes against analytic solutions. Then, we show the application of the numerical schemes to the integration of test particle trajectories in general relativistic magnetohydrodynamic (GRMHD) simulations, by modifying the algorithms to handle grid-based electromagnetic fields. We test this approach by simulating ensembles of charged particles in a static GRMHD configuration obtained with the Black Hole Accretion Code (\texttt{BHAC}).