The ensemble density functional theory is valuable for simulations of metallic systems due to the absence of a gap in the spectrum of the Hamiltonian matrices. Although the widely used self-consistent field iteration method can be extended to solve the minimization of the total energy functional with respect to orthogonality constraints, there is no theoretical guarantee on the convergence of these algorithms. In this paper, we consider an equivalent model with a single variable and a single spherical constraint by eliminating the dependence on the fractional occupancies. A proximal gradient method is developed by keeping the entropy term but linearizing all other terms in the total energy functional. Convergence to the stationary point is established. Numerical results in the KSSOLV toolbox under the Matlab environment show that they can outperform SCF consistently on many metallic systems.