The material point method (MPM) simulation of large-scale landslides can provide useful information to response to landslide hazards in mountainous regions. However, the structured background meshes commonly adopted in MPM fail to accurately capture the complex terrain above which the sliding mass moves. In this work, a novel terrain contact algorithm is proposed and integrated into a GPU-accelerated 3D MPM to model large-scale landslides. The complex terrain is modelled by a fixed triangular surface mesh generated according to the digital elevation model (DEM) and the interaction between the sliding mass and the basal surface is modelled by a contact algorithm between the triangular surface mesh and the material points. Furthermore, a contact detection algorithm based on level-set functions is developed to improve the contact detection efficiency. The detailed computational formulations and numerical implementation based on GPU parallel computing are presented. After the correctness and the efficiency being verified by two numerical examples, the proposed method is applied to simulate the catastrophic Baige landslide that occurred in China in 2018. The results indicate that the present method is both accurate and efficient and can provide a useful tool for numerical study on landslides.