A new procedure for the optimization of the exponents, alpha(j), of Gaussian basis functions, Y-l(m)(theta,phi)r(l)e(j)(-alpha)r(2), is proposed and evaluated. The direct optimization of the exponents is hindered by the very strong coupling between these nonlinear variational parameters. However, expansion of the logarithms of the exponents in the orthonormal Legendre polynomials, P-k, of the index, j: ln alpha(j)=Sigma(k=0)(kmax)A(k)P(k)((2j-2)/(N-prim-1)-1), yields a new set of well-conditioned parameters, A(k), and a complete sequence of well-conditioned exponent optimizations proceeding from the even-tempered basis set (k(max)=1) to a fully optimized basis set (k(max)=N-prim-1). The error relative to the exact numerical self-consistent field limit for a six-term expansion is consistently no more than 25% larger than the error for the completely optimized basis set. Thus, there is no need to optimize more than six well-conditioned variational parameters, even for the largest sets of Gaussian primitives. (C) 2003 American Institute of Physics.