Elements heavier than iron are being created in massive, explosive astrophysical scenarios such as core-collapse supernovae and neutron star mergers. Heavy element nucleosynthesis is a very complex process that requires the knowledge of the properties, mainly nuclear masses, neutron capture and beta-decay rates, of thousands of nuclei. However, due to the limitations of current experimental facilities, only a relatively small number of nuclei have so far been studied experimentally. Given the amount of nuclei, and the regions of the nuclear chart involved in the heavy element nucleosynthesis, it is necessary to use models that can be reliably applied to even the most exotic nuclei. Thus, the use of microscopic nuclear structure models presents itself as the logical choice, because the underlying microscopic theory of nuclear interaction is expected to be valid across the whole nuclear chart. In this manuscript, we present the results of a large-scale calculation of β-decay rates using a fully self-consistent theoretical framework based on the relativistic nuclear energy density functional. Taking into account the first-forbidden transitions, we are able to determine the regions of the nuclear chart where these transitions are critical for the description of the decay rates. Finally, we examine the β-delayed neutron emission in all neutron-rich nuclei, from the valley of stability to the limits of nuclear binding. The results of the study are applied in an r-process simulations where we obtain an improved description of the r-process abundance pattern.