Should we first relax the structure (Volume, cell shape and ions) with non magnetic ground state and then use ISPIN=2 for single point geometry calculation with the relaxed structure? Would this method give more stable ground state energy?
It depends on the specifics of your system. If you're trying to do this for a compound with small local magnetic moments, then it may be a clever choice to first do the optimization using a non-magnetic setting and then complete it with a a fully spin-polarized one (ISPIN=2). On the other hand, if you're working on strongly correlated materials containing high spin-multiplicity elements such as 3d transition metals or rare-earth f-electron elements such as Sm, Eu or Gd, you may never get a converged solution using ISPIN=1. The reason is, such localized states want to obey Hund's exclusion rule and get spin-polarized which obviously can't be fulfilled with this setting. As a result, they keep fluctuating around the Fermi-level, leading to instabilities or even divergence in the SCF loop.