In drug discovery, molecular dynamics (MD) simulation for protein-ligand binding provides a powerful tool for predicting binding affinities, estimating transport properties, and exploring pocket sites. While significant strides have been made in advancing more efficient MD simulations, the accurate modeling of extended-timescales simulations remains a considerable challenge. To address this issue, we propose NeuralMD, the first approach to learning to predict binding dynamics. It builds upon a novel multi-grained group symmetric framework and effectively incorporates the physics laws from two perspectives: (1) The geometric representation of the protein-ligand complex should be SE(3)-equivariant and can sufficiently capture the particle interplay between protein and ligand, and (2) The trajectory learning for binding dynamics needs to follow Newtonian mechanics. To verify the effectiveness of NeuralMD, we design ten single-trajectory and three multi-trajectory binding simulation tasks. Quantitatively, NeuralMD outperforms three competitive machine learning baselines for binding MD simulation.