Based on a recently proposed exponential parametrization of the one-electron atomic-orbital density matrix of a single-determinant wave function [Chem. Phys. Lett. 327, 397 (2000)], we present an implementation of the direct optimization of the atomic-orbital density matrix as an alternative to the diagonalization of the Fock-Kohn-Sham matrix when solving the Roothaan-Hall self-consistent field equations. The optimization of the density matrix is carried out by the conjugate-gradient method with a multilevel nondiagonal preconditioner and is well suited to linear scaling. Although a diagonal preconditioner may be sufficient for minimal basis sets and large highest occupied molecular orbital-lowest unoccupied molecular orbital (HOMO-LUMO) gaps, a nondiagonal preconditioner is needed in more difficult cases-that is, for basis sets with polarization and diffuse functions and for systems with small HOMO-LUMO gaps. Redundancies of the exponential parametrization of the density matrix are handled by a projection technique, thereby avoiding singular equations in the optimization of the density matrix.