We develop and analyze Riemannian optimization methods for computing ground states of rotating multicomponent Bose?Einstein condensates, modeled as minimizers of the Gross?Pitaevskii energy functional. To address the non-uniqueness of ground states induced by phase invariance, we work on a quotient manifold that factors out the symmetry and endow it with a general Riemannian metric tailored to the problem?s energy landscape. Using an auxiliary phase-aligned iteration and fixed-point convergence theory, we establish a general local convergence framework for Riemannian gradient descent methods with explicit contraction rates. Specializing to two metrics, we derive the energy-adaptive method (eaRGD), which ensures monotone energy decay and global convergence, and the Lagrangian-based method (LagrRGD), for which we provide the first rigorous local convergence proof for Gross?Pitaevskii energy functionals, valid in both rotating and non-rotating settings. Incorporating second-order information, LagrRGD achieves faster local convergence than eaRGD. Numerical experiments confirm the theoretical results and demonstrate the efficiency of both approaches.