The study of ion transport in electrochemically active materials for energy storage systems requires simulations on quantum- atomistic- and mesoscales. The methods accessing these scales not only have to be effective but also well compatible to provide a full description of the underlying processes. We propose to adapt the Generalized Shadow Hybrid Monte Carlo (GSHMC) method to atomistic simulation of ion intercalation electrode materials for batteries. The method has never been applied to simulations in solid state chemistry but it has been successfully used for simulation of biological macromolecules, demonstrating better performance and accuracy than can be achieved with the popular molecular dynamics (MD) method. It has been also extended to simulations on meso-scales, making it even more attractive for simulation of battery materials. We combine GSHMC with the dynamical Core-Shell model to incorporate polarizability into the simulation and apply the new Modified Adaptive Integration Approach, MAIA, which allows for a larger time step due to its excellent conservation properties. Also, we modify the GSHMC method, without losing its performance and accuracy, to reduce the negative effect of introducing a shell mass within a dynamical shell model. The proposed approach has been tested on olivine NaFePO4, which is a promising cathode material for Na-ion batteries. The calculated Na-ion diffusion and structural properties have been compared with the available experimental data and with the results obtained using MD and the original GSHMC method. Based on these tests, we claim that the new technique is advantageous over MD and the conventional GSHMC and can be recommended for studies of other solid-state electrode and electrolyte materials whenever high accuracy and efficient sampling are critical for obtaining tractable simulation results.