Simulating the deformation behavior of crystals provides fundamental insights into the mechanics of polycrystalline materials such as metals and alloys. Single-crystal plasticity models, which are based on the crystallographic structure of individual grains, describe deformation at the continuum level. By invoking the principle of maximal dissipation, the problem can be reformulated as a nonlinear, non-convex optimization problem. For rate-independent formulations, algorithmic challenges arise from the non-uniqueness of active slip systems and the potential for multiple solutions, necessitating robust and efficient computational methods. Conventional approaches include regularization via the visco-plastic limit case, simplifications to ensure uniqueness, or active-set search algorithms to identify active slip systems. In recent years, new methods have emerged to address this constrained optimization problem. Among these, the interior point method—which penalizes constraint violations using logarithmic barrier functions—has shown significant promise. In this work, we present a robust and efficient interior point formulation for finite strain crystal plasticity [1, 2], specifically designed to handle complex hardening scenarios, including highly nonlinear hardening functions and self-/latent hardening mechanisms. The performance of the algorithm is evaluated through numerical examples, demonstrating its reliability and computational efficiency.