The determination of accurate nuclear matrix elements for ββ decay processes is a challenge for nuclear theory and can have a strong impact in neutrino physics. Large Scale Shell Model (LSSM) calculations are among the best tools for such determination and recent developments have allowed to extend its application domains. In particular, systematic studies of nuclear matrix elements calculations have been now undertaken in this framework for most of the ββ emitters. These calculations are crucial in the determination of the most favorable emitters in the forthcoming generation of ββ experiments. The present paper focuses on the recent advances and remaining difficulties of shell model calculations for the neutrinoless mode. Stability and predictive power of the results will be discussed.