Molecular beam epitaxy (MBE) is an important and challenging research topic in material science. In this paper, we propose a new fully discrete scheme for the well-celebrated continuum MBE model with slope selection. First of all, we use a multi-step strategy to discretize the MBE model in time. The obtained semi-discrete scheme is proved to possess properties of total mass conservation, unconditionally energy stability and uniquely solvability. The rigorous error estimate is then conducted to show its second-order convergence. The semi-discrete scheme is further discretized in space using the Fourier pseudo-spectral method. The fully discrete scheme is also shown to preserve mass-conservation and energy-dissipation properties. Afterward, several numerical examples are presented to validate the accuracy and efficiency of our proposed scheme. In particular, the scaling law for the roughness growing and effective energy decaying are captured during long-time coarsening dynamic simulations. The idea proposed in this paper could be readily utilized to design accurate and stable numerical approximations for many other energy-based phase field models.