The new techniques in adaptive optics, free space optical(FSO) communication rely on the use of numerical simulations for atmospheric turbulence to evaluate the performance of the system. The simulation of turbulence phase screen is the heart of numerical simulations which produces random wavefront phase perturbations with the correct statistical properties corresponding to models of optical propagation through atmospheric turbulence. The phase-screen simulation techniques can be roughly divided into fast Fourier transform (FFT) method and matrix-based method. Because of a better performance in computation time, the FFT method is generally used for modeling the performance of a real system. But the classical FFT method has a main deficiency of oversample in low frequency region, which leads to the lost of accuracy. To overcome this deficiency, many methods have been proposed for compensating for the oversample of low frequency components, in the last decades. Essentially, these methods achieve a higher accuracy at the expense of computation time. A good compensation method should take into consideration both accuracy and computation time. 　　To achieve higher accurcy and lower computational cost simultaneously, we develop a hybrid method to generate turbulence phase screen, i.e. the classical FFT model is mixed with the sparse spectrum model. We first extract the low frequency region from the frequency grid of FFT model, and resample this region with 16 samples. It is found that the accuracy of phase screen is related to the distribution of these samples, and there must be an optimum distribution that can minimize the relative error between expected structure function and theoretical structure function in the low frequency region. So it permits one to use optimization algorithm to find the optimized distribution of low frequency samples. Here an improved gravity search algorithm is adopted in which the memory of each particle is taken into consideration. The optimization parameters are determined after a lot of tests, and the robustness testing shows that the algorithm is effective. To compare with existing subharmonic method, we choose the same parameters of phase screen as those used in the expanded subharmonic method, generate 1000 phase screens for each method, compute the phase structure function, and we also compare our results with those from the theoretical structure function. The comparison result shows that the curve of phase structure function generated by our method is nearly consistent with the theoretical one, the maximum relative error in low frequency region is about 0.063% which is much better than that from the expanded subharmonic method 5%. Finally in this paper, the computational cost is analyzed, showing that the generation speed for our method is at least 4.5 times as fast as that for the Johansson’s method.