An efficient method for optimizing single-determinant wave functions of medium and large systems is presented. It is based on a minimization of the energy functional using a new set of variables to perform orbital transformations. With this method convergence of the wave function is guaranteed. Preconditioners with different computational cost and efficiency have been constructed. Depending on the preconditioner, the method needs a number of iterations that is very similar to the established diagonalization-DIIS approach, in cases where the latter converges well. Diagonalization of the Kohn-Sham matrix can be avoided and the sparsity of the overlap and Kohn-Sham matrix can be exploited. If sparsity is taken into account, the method scales as O(MN2), where M is the total number of basis functions and N is the number of occupied orbitals. The relative performance of the method is optimal for large systems that are described with high quality basis sets, and for which the density matrices are not yet sparse. We present a benchmark calculation on a DNA crystal containing 2x12 base pairs, solvent and counter ions (2388 atoms), using a TZV(2d,2p) basis (38 688 basis functions) and conclude that the electronic structure of systems of this size can now be studied routinely.