In this paper, we develop an iterative scheme to construct multiscale basis functions within the framework of the Constraint Energy Minimizing Generalized Multiscale Finite Element Method (CEM-GMsFEM) for the mixed formulation. The iterative procedure starts with the construction of an energy minimizing snapshot space that can be used for approximating the solution of the model problem. A spectral decomposition is then performed on the snapshot space to form global multiscale space. Under this setting, each global multiscale basis function can be split into a non-decaying and a decaying parts. The non-decaying part of a global basis is localized and it is fixed during the iteration. Then, one can approximate the decaying part via a modified Richardson scheme with an appropriately defined preconditioner. Using this set of iterative-based multiscale basis functions, first-order convergence with respect to the coarse mesh size can be shown if sufficiently many times of iterations with regularization parameter being in an appropriate range are performed. Numerical results are presented to illustrate the effectiveness and efficiency of the proposed computational multiscale method.