The Moving Particle Semi-implicit (MPS) method performs well in simulating violent free surface flow and hence becomes popular in the area of fluid flow simulation. However, the implementations of searching neighbouring particles and solving the large sparse matrix equations (Poisson-type equation) are very time-consuming. In order to utilize the tremendous power of parallel computation of Graphics Processing Units (GPU), this study has developed a GPU-based MPS model employing the Compute Unified Device Architecture (CUDA) on NVIDIA GTX 280. The efficient neighbourhood particle searching is done through an indirect method and the Poisson-type pressure equation is solved by the Bi-Conjugate Gradient (BiCG) method. Four different optimization levels for the present general parallel GPU-based MPS model are demonstrated. In addition, the elaborate optimization of GPU code is also discussed. A benchmark problem of dam-breaking flow is simulated using both codes of the present GPU-based MPS and the original CPU-based MPS. The comparisons between them show that the GPU-based MPS model outperforms 26 times the traditional CPU model. moving particle semi-implicit method (MPS), graphics processing units (GPU), compute unified device architecture (CUDA), neighbouring particle searching, free surface flow PACS: 47.10.ad, 45.10Db In the community of Computational Fluid Dynamics, many mesh-based methods have been developed so far to solve various problems of fluid flows. However, the well developed numerical methods such as the finite difference (volume) method, finite element method and boundary element method remain great challenges in dealing with the violent free surface flow with large deformation. Meshless methods were developed to overcome the inherent difficulties in the numerical methods based on mesh generation. Firstly, Monaghan [1] proposed the Smoothed Particle Hydrodynamics (SPH) method to deal with the astrophysical hydrodynamic problems. In SPH, the kernel approximations are used to interpolate the unknowns using particles in the Lagrangian sense. Later, this method was extended to computational fluid dynamics. Monaghan [2] and Shao [3] improved the original SPH to simulate the incompressible Newtonian fluid flows with a free surface. In addition, Koshizuka [4] proposed another particle-based method called Moving Particle Semi-implicit (MPS) for the incompressible viscous fluid flow. By MPS, the fluid flow is represented by moving particles and the convection term is calculated by the motion of these moving particles. Even if the fragmentation or merging of fluid occurs, the interfaces can be represented clearly. MPS method has been used successfully in solving many problems, such as the breaking waves [5], water sloshing [6] and tracking the free surface [7]. However, just as other particle-based CFD methods, MPS is shown to be very time-consuming, mainly because of the required searching procedure for the neighbouring