Contradicting conclusions have been reported regarding the change of Li diffusivity as a function of Li concentration. For example, using the galvanostatic intermittent titration Technique (GITT) Park et al. found a phase transition dependent Li diffusivity versus Li content and an overall increasing Li diffusivity with a higher Li content in Li x C 6 while Piao et al. showed a monotonic relation between Li diffusivity versus Li content and a decreasing Li diffusivity with a higher Li content using AC impedance measurements. Computational simulations have also been implemented to calculate Li diffusion coefficients in graphite. Using first-principles calculations combined with Monte Carlo simulations, Persson et al. reported a high Li diffusivity of 10 -7 to 10 -6 cm 2 /s in the direction parallel to the graphene layers. In their calculations, the interlayer spacing between empty graphene layers was fixed during simulations while experimentally, the layer spacing increases with a higher Li content. In addition, only stage 1 and 2 structures are considered and Li diffusivity was calculated for x > 0.2 in Li x C 6 . In this study, we investigate the phase transitions and Li diffusion in the Li x C 6 system as a function of stage structure, Li content, and temperature. Given the high computational cost of ab initio molecular dynamics due to the large size and time scales required to obtain sufficient diffusion statistics, especially at low temperatures, we train a machine learning force field (MLFF) based on our recently reported charge recursive neural network architecture (QRNN) for simulating the Li diffusion in the Li x C 6 system with 0 ≤ x ≤ 1 and stage 1, 2, 3, 4, and 1L. Our QRNN model was trained with a dataset labeled by density functional theory (DFT) calculations with the nonlocal van der Waals functional rVV10. We demonstrate that comparing with the PBE functionals and previously reported ab initio results, the rVV10 functional accurately reproduces the experimental lattice parameters, phase stability, and voltage profile of the Li x C 6 system. Previously, we developed a QRNN model trained to non-periodic DFT data for liquid electrolyte simulations and our predicted bulk thermodynamic and transport properties of pure carbonate solvents and electrolyte mixtures agree well with experimental results. In this work, we demonstrate that our QRNN architecture also works well for periodic systems. Our QRNN model is accurate and stable, allowing us to perform large scale molecule dynamics simulations to explore Li intercalation and diffusion in the graphite anode.