Conventional semiempirical methods using diagonalization are not practical for calculations on molecular systems-containing more than a feu hundred atoms because of O(N-3) time and O(N-2) memory requirements, where N is the number of atoms. Currently, the time dominating step is diagonalization of the Fock matrix. This paper demonstrates how O(N-3) diagonalization and O(N-2) memory requirements are eliminated by using a conjugate gradient search for the density matrix with sparse matrix techniques. Our method makes high accuracy energy calculations on molecules containing thousands of atoms possible on the typical workstation. Benchmark examples are presented on polyglycine chains (20000 atoms), water clusters (up to 1800 atoms), and nucleic acids (up to 6304 atoms).