In geotechnical applications the success of the discrete element method (DEM) in simulating fundamental aspects of soil behaviour has increased the interest in applications for direct simulation of engineering scale boundary value problems (BVP’s). The main problem is that the method remains relatively expensive in terms of computational cost. A non-negligible part of that cost is related to specimen creation and initialization. As the response of soil is strongly dependant on its initial state (stress and porosity), attaining a specified initial state is a crucial part of a DEM model. Different procedures for controlled sample generation are available. However, applying the existing REV-oriented initialization procedures to such models is inefficient in terms of computational cost and challenging in terms of sample homogeneity. In this work a simple but efficient procedure to initialize large-scale DEM models is presented. Periodic cells are first generated with a sufficient number of particles matching a desired particle size distribution (PSD). The cells are then equilibrated at low-level isotropic stress at target porosity. Once the cell is in equilibrium, it is replicated in space in order to fill the model domain. After the domain is thus filled a few mechanical cycles are needed to re-equilibrate the large domain. The result is a large, homogeneous sample, equilibrated under prescribed stress at the desired porosity. The method is applicable to both isotropic and anisotropic initial stress states, with stress magnitude varying in space.