An interpolation method is presented for calculating a divergence-free magnetic field at arbitrary locations in space from a representation of that field on a discrete grid. This interpolation method is used along with symplectic integration to perform particle trajectory integrations with good conservation properties. These integrations are better at conserving constants of motion and adiabatic invariants than standard, nonsymplectic Runge-Kutta integration schemes. In particular, we verify that carrying out particle integrations with interpolated magnetic fields that satisfy $nabla$ ¿ B = 0 yields a better conservation of the first adiabatic invariant ¿. Comparisons are made between the different integration methods for proton trajectories in an ideal dipole magnetic field and in the cusp region of the magnetosphere as determined numerically from a global magnetohydrodynamic (MHD) model for realistic solar wind conditions. In the case considered, we find that particle trapping can occur in the cusp only if the self-consistent electric field obtained from the MHD code is not taken into account. When that electric field is included in the calculation of particle trajectories, we fail to find any particle trapping in the cusp region.