We present generalized loading-unloading contact laws for elasto-plastic spheres with bonding strength. The proposed mechanistic contact laws are continuous at the onset of unloading by means of a regularization term, in the spirit of a cohesive zone model, that introduces a small and controllable error in the conditions for interparticle breakage. This continuity property is in sharp contrast with the behavior of standard mechanistic loading and unloading contact theories, which exhibit a discontinuity at the onset of unloading when particles form solid bridges during plastic deformation. The formulation depends on five material properties, namely two elastic properties (Young's modulus and Poisson's ratio), two plastic properties (a plastic stiffness and a power-law hardening exponent) and one fracture mechanics property (fracture toughness), and its predictions are in agreement with detailed finite-element simulations. The numerical robustness and efficiency of the proposed formulation are borne out by performing three-dimensional particle mechanics static calculations of microstructure evolution during the three most important steps of powder die-compaction, namely during compaction, unloading, and ejection. These simulations reveal the evolution, up to relative densities close to one, of microstructural features, process variables and compact mechanical attributes which are quantitatively similar to those experimentally observed and in remarkable agreement with the (semi-)empirical formulae reported in the literature.