A numerical framework is proposed and analyzed for computing the ground state of Bose-Einstein condensates. A gradient flow approach is developed, incorporating both a Lagrange multiplier to enforce the L^2 conservation and a free energy dissipation. An explicit approximation is applied to the chemical potential, combined with an exponential time differencing (ETD) operator to the diffusion part, as well a stabilizing operator, to obtain an intermediate numerical profile. Afterward, an L^2 normalization is applied at the next numerical stage. A theoretical analysis reveals a free energy dissipation under a maximum norm bound assumption for the numerical solution, and such a maximum norm bound could be recovered by a careful convergence analysis and error estimate. In turn, the proposed method preserves the following theoretical properties: (1) an explicit computation at each time step, (2) unconditional free energy dissipation, (3) L^2 norm conservation at each time step, (4) a theoretical justification of convergence analysis and error estimate. Comprehensive numerical experiments validate these theoretical results, demonstrating excellent agreement with established reference solutions.
