A theory of globally convergent trust-region methods for self-consistent field electronic structure calculations that use the density matrices as variables is developed. The optimization is performed by means of sequential global minimizations of a quadratic model of the true energy. The global minimization of this quadratic model, subject to the idempotency of the density matrix and the rank constraint, coincides with the fixed-point iteration. We prove that the global minimization of this quadratic model subject to the restrictions and smaller trust regions corresponds to the solution of level-shifted equations. The precise implementation of algorithms leading to global convergence is stated and a proof of global convergence is provided. Numerical experiments confirm theoretical predictions and practical convergence is obtained for difficult cases, even if their geometries are highly distorted. The reduction of the trust region is performed by a strategy that uses the structure of the energy function providing the algorithm with a nice practical behavior. This framework may be applied to any problem with idempotency constraints and for which the derivative of the objective function is a symmetric matrix. Therefore, application to calculations based both on Hartree-Fock or Kohn-Sham density functional theory are straightforward.