ABSTRACT We present a new magnetohydrodynamic-particle-in-cell (MHD-PIC) code integrated into the athena++ framework. It treats energetic particles as in conventional PIC codes, while the rest of thermal plasmas are treated as background fluid described by MHD, thus primarily targeting at multiscale astrophysical problems involving the kinetic physics of the cosmic rays (CRs). The code is optimized towards efficient vectorization in interpolation and particle deposits, with excellent parallel scaling. The code is also compatible with static/adaptive mesh refinement, with dynamic load balancing to further enhance multiscale simulations. In addition, we have implemented a compressing/expanding box framework that allows adiabatic driving of CR pressure anisotropy, as well as the δf method that can dramatically reduce Poisson noise in problems where distribution function f is only expected to slightly deviate from the background. The code performance is demonstrated over a series of benchmark test problems, including particle acceleration in non-relativistic parallel shocks. In particular, we reproduce the linear growth of the CR gyro-resonant (streaming and pressure anisotropy) instabilities, under both the periodic and expanding/compressing box settings. We anticipate the code to open up the avenue for a wide range of astrophysical and plasma physics applications.