Metachronal motions are ubiquitous in terrestrial and aquatic organisms and have attracted substantial attention in engineering for their potential applications. Hard-magnetic soft materials are shown to provide new opportunities for metachronal wave-modulated robotic locomotion by multi-agent active morphing in response to external magnetic fields. However, the design and optimization of such magnetic soft robots can be complex, and the fabrication and magnetization processes are often delicate and time-consuming. Herein, a computational model is developed that integrates granular models into a magnetic-lattice model, both of which are implemented in the highly efficient parallel computing platform large-scale atomic/molecular massively parallel simulator (LAMMPS). The simulations accurately reproduce the deformation of single cilium, the metachronal wave motion of multiple cilia, and the crawling and rolling locomotion of magnetic cilia soft robots. Furthermore, the simulations provide insight into the spatial and temporal variation of friction forces and trajectories of cilia tips. The results contribute to the understanding of metachronal wave-modulated locomotion and potential applications in the field of soft robotics and biomimetic engineering. The developed model also provides a versatile computational framework for simulating the movement of magnetic soft robots in realistic environments and has the potential to guide the design, optimization, and customization of these systems.